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ABSTRACT 

A  method  is  developed  in  terms  of  the  Cauer  Second  Form 
representation  of  continued  fractions  as  a  means  of  designing 
linear  single-input  single-output  (SISO)  control  systems. 
Optimal  closed  loop  solutions  corresponding  to  a  set  of 
prescribed  eigenvalues  are  obtained  through  minimization  of 
a  quadratic  performance  index.   The  partitioning  method  of 
the  Cauer  Second  Form  for  system  simplification  is  presented 
with  a  simplified  inversion  technique  for  the  reduced  order 
system. 
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I.   INTRODUCTION 

The  purpose  of  this  research  was  to  develop  an  algorithm 
for  obtaining  optimal  closed  loop  solutions  corresponding 
to  a  set  of  prescribed  eigenvalues  for  single-input  single- 
output  (SISO)  control  systems.   It  was  desired  that  the 
algorithm  be  adaptable  to  digital  computer  techniques  and 
unrestricted  by  system  order. 

The  Cauer  Second  Form  for  system  dynamics  representation 
was  chosen  over  other  alternatives  because  of  the  regular 
pattern  of  the  state  and  output  matrices,  and  the  method 
of  linear  system  simplification. 

In  Chapter  II,  several  basic  properties  of  both  Cauer 
First  and  Second  Forms  are  presented  from  the  theory  of 
continued  fractions.   A  simple  and  efficient  algorithm  is 
also  developed  for  inversion  of  the  continued  fraction  in 
either  form,  independent  of  Routh ' s  algorithm. 

In  Chapter  III,  the  method  of  linear  system  order 
reduction  based  on  the  Cauer  Second  Form  is  amplified. 
The  emphasis  on  this  area  was  primarily  to  elucidate  the 
various  methods  previously  employed  for.  system  simplifi- 
cation. 

The  original  theoretical  work  of  this  thesis  is  pre- 
sented in  Chapter  IV.   The  objective  was  to  obtain  closed 
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loop  solutions  corresponding  to  a  prescribed  set  of  eigen- 
values.  While  minimizing  a  certain  cost  function,  which 
met  desired  system  characteristics.   It  is  shown,  by  examples, 
that  the  derived  algorithm  is  equally  capable  of  handling 
systems  with  multiple  and/or  complex,  as  well  as,  unique 
sets  of  real  eigenvalues. 

The  final  chapter,  Chapter  V,  presents  a  discussion  of 
results  and  suggests  areas  for  future  study. 


II.   PROPERTIES  OF  CAUER  FIRST  AND  SECOND  FORMS 


A.   CLOSED  LOOP  SYSTEM  IN  CAUER  FIRST  AND  CAUER  SECOND 
FORMS 

Consider  the  closed  loop  transfer  function  given  by 

n-1 

I   b.s1 

Y(s)     L=0   x 

U(s)        n-1      .  ' 
n  ,   r      1 
s   +   Z   a .  s 

i=0  "■ 


(2-1) 


with  block  diagram  as  given  in  Figure  2.1.   Equation  (2-1) 
can  be  expanded  into  the  Cauer  Forms  of  continued  fractions 
as  follows. 

1.   Cauer  First  Form 

a.  Arrange  the  numerator  and  denominator  poly- 
nomials in  descending  order. 

b.  Perform  continued  division. 


■u  n-1  ,  „n- 2  .  .  -u  -L 
w  x  b  ,s  +  b  0s  +  ...  +  b,  s  +  b 
Y(  s)  _   n-1 n-2 1 o 

U(s)    n  .  a    on-l  .  a     on-2  .      .  a   .a 
s   +a  -,    s     +a   ns     +...+a,  s  +  a 
n-1         n-2  1   o 


(2-2) 


hxs  + 


h0  + 1 

2    * (2-3) 


h3s  + 


\    + 


10 
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2  .   Cauer  Second  Form 

a.   Invert  the  numerator  and  denominator  and  arrange 
the  polynomials  in  ascending  order. 


n-2         n-1     n 
„,  v    a   +  a.  s  +  ....a  ns           +  a   ,s     +  s 
Y(s)  _  _o 1 n-2 n-1 

U(s)    .     ,  „  .  .,     n-2    .     n-1 

bo  +  blS  +  +bn-2S     +  bN-lS 


b.   Perform  continued  division. 
Y(s)  _ 


(2-4) 


UTsT 


hl  + 


h2 

+  1  (2-5) 


3 


h3  + 


\   + 


or 


hi  + 


(2-6) 


h  + 
H2 


h3  + 


\    + 


Block  diagrams  of  both  systems  are  shown  in  Figures  2.2  and 
2.3. 
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B.   PHYSICAL  INTERPRETATION  OF  DOMINANT  TERMS 

RESULTING  FROM  CONTINUED  FRACTION  EXPANSION 

It  is  known  that  the  most  dominant  terms  in  equations 
(2-3)  and  (2-5)  are  the  first  quotients,  h. s  and  h, ,  re- 
spectively.  A  meaningful  interpretation  for  these  terms 
can  be  accomplished  by  applying  the  initial  value  and  final 
value  theorems.    Letting  Y(s)/U(s)  =  F(s),  by  an  asymptotic 
expansion  approximation: 

1.   For  Cauer  First  Form 

lim  f(t)  ~  lim  sF(s)  s  ~  (2-7) 

hl 
t->o        s^-°° 


and 


lim  f(t)  ~  lim  sF(s)  ~  h2  +  h^.  (2-8) 


+ 


2.   For  Cauer  Second  Form 


lim  f(t)  ~  lim  sF(s)  a  h2  +  h^  (2-9) 


lim  f(t)  =  lim  sF(s)  =  h,.  (2-10) 

t-*»        s-*o 


lim  f(t)  must  exist. 
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Equations  (2-7)  and  (2-10)  are  of  considerable  interest  since 
they  involve  the  dominant  term,  h, .   The  implication  is  that 
the  Cauer  First  Form  emphasizes  the  initial  or  transient 
response  of  the  system;  whereas,  the  Cauer  Second  Form  em- 
phasizes the  final  or  steady  state  response  of  the  system. 
In  general,  the  quotients  lower  in  position  in  the  continued 
fraction  expansion  have  less  influence  on  the  performance  of 
the  system  as  a  whole,  (h.  has  less  influence  than  h.,  where 
i<j).   Because  many  systems  must  meet  a  set  of  steady  state 
conditions,  the  Cauer  Second  Form  will  be  used  for  the 
prescribed  eigenvalue  problem. 

C.   CONTINUED  FRACTION  INVERSION 

The  theory  of  continued  fractions  was  first  associated 
with  Routh's  Algorithm  by  Wall  in  1945,  [1]  and  [2].   The 
following  year  Frank  [3]  extended  and  modified  Wall's  work 
to  include  complex  coefficients.   Both,  however,  applied 
Routh's  algorithm  only  to  continued  fraction  expansions, 
not  to  the  problem  of  inversion. 

In  1969,  Chen  and  Shieh  [4]  developed  an  algorithm 
method  for  converting  a  continued  fraction  into  a  rational 
fraction  of  two  polynomials.   Their  method,  which  makes  use 
of  Routh's  algorithm,  is  presented  below. 

If  the  elements,  h.,  are  known  for  any  continued 
fraction,  then  the  state  and  output  equations  can  be 
written  immediately  from  Figures  2.2  or  2.3. 
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n 


h2hi        huhi 


h6hl 


2n    1 


h0h,        h.-CliT+ha)        hc(h   +h0)         ...h0    (h   +h.) 
21  413  613  2nl3 


h2hx        VW         h6(hl+h3+h5) 


h2hx        hlf(h1+h3)        h6(h1+h3+h5)...h2n(h1+..+h2n.1) 


z. 

1 

1 

Z2 

1 

Z3 

+ 

1 

z 

1 

n  _ 

(2-11) 


z    =    Hz    +    Dr 


(2-12) 


C    =    [h2   \   h6    •'•    h2x] 


n 


(2-13) 


C    =   L   z 


(2-14) 


17 


The  coefficients  of  the  characteristic  polynomial, 
|sI-H|,  become  the  first  row  elements  of  the  required  Routh 
array.   The  next  sequence  of  steps  in  determining  the  (j+2)th 
row  of  the  Routh  array  is  to  successively  let 


h.  =  0 
3 


and 


hj+l  =  ° 


(2-15) 


for  js[l,3,5, . . .2n_l] 

and  evaluating  the  remaining  (n-k)x£ n-k)  "ft."  matrix,  where 
k  =  (j+l)/2,  up  to  k=n-l,  i.e.,  for  n  arbitrary  and  j=l; 
the  3rd  row  of  the  Routh  array  becomes  (after  h,  and  h„  are 
set  equal  to  zero)  the  coefficients  of 


sI-Hj, 


where  the  (n-l)x(n-l)  H,  matrix  is 


h.  ,h_      h_h_ 
4    3         6    3 


h4h3      VW 


h2nh3 


h„h0      hc(h0+hc)    h0    (h   +hc) 

M-3635  2n3o 


,h0^(h0+  .  .  -h0^    -,  ) 
2n      3  2n-l 


(2-16) 


This  process  repeats  until  the  system  state  matrix  is 
reduced  to  a  single  element,  H„   .,,  yielding  the  (2n-l)th 
row  in  Routh ' s  array.   It  is  observed  that  each  successive 
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odd  numbered  row  contains  one  less  element  than  it's 
predecessor.   By  inserting  leading  zeros  in  the  3rd,  5th, 
...,  (2n+l)th  row,  the  matrix,  P,  is  formed. 

3rd 
5  th 
7th       0       0      P33   1  (2-17) 


(2n+l)th 


Pll 

P12 

P13   

1 

0 

P22 

P23   

1 

0 

0 

P33   

1 

0 


0 


0 


The  matrix,  P,  is  the  linear  transformation  matrix  required 
to  obtain  a  linear  system  in  Cauer  Second  Form  from  phase 
variable  (canonical)  form.   Continuing,  the  second  row  of 
the  Routh  array  is  obtained  from  the  output  matrix,  L, 
and  the  above  transformation: 

c  =  Lz    (continued  fractions) 

z  =  Px    (linear  transformation) 

y  =  Cx    (output  equation,  phase  variable  form) 

Therefore , 


C  =  L  P  . 


(2-18) 


C  is  an  (lxn)  vector  whose  elements  are  the  seoond  row  of 
the  Routh  array. 

Consider  the  Routh  array  as  an  (n+l)x(2n+l)  matrix  with 
typical  element  r...   The  quotients,  h.,  of  the  continued 
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fraction  expansion  can  be  expressed  as: 

h±    =  — ii-  (2-19) 

1   ri+l,l 

From  this  relationship  and  knowledge  of  how  the  Routh 

array  is  generated,  the  remaining  even  numbered  rows  of  the 

array  can  be  found.   The  transfer  function  as  a  ratio  of  two 
polynomials  is  written  as: 


n        .  , 

T(s)  =  ^~ (2-20) 

n^        i-1 
Z   r,  .  s 

L=l   1»3 


Chen  and  Shieh  [4]  contend  that  this  method  is  the 
easiest  in  attaining  the  inversion.   The  author  disagrees 
and  presents  a  simpler  iterative  method  based  on  the  in- 
version technique  for  the  Generalized  Cauer  Form  given  by 
Goldman  [5].   The  method  is  equally  suited  to  both  Cauer 
First  and  Cauer  Second  Forms,  requiring  no  prior  knowledge 
of  Routh's  algorithm.   Assuming  all  h.'s  are  known,  and 
non-zero,  in  equation  (2-3)  or  (2-5),  let: 


a.  =  h2i_x  (2-21) 

b.  =  h„.  (2-22) 
l     2i 

for  i  e.[l ,  2  ,  .  .  .  .  ,n]  . 
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1 .   Inversion  of  Cauer  First  Form 

Initialize  two  (n+lxl)  vectors  C  and  D 


C  =  [cQ  c1  c2  cr]  (2-23) 

D  =  [dQ  d1  d2  dR]  (2-24) 


to  all  zeros ,  except 


cn  =  bn  (2-25) 

dn-l  =  ax  X  Cn  (2"26) 

dn  =  1.  (2-27) 


The  following  set  of  equations  are  first  solved  for  i=l. 


c      ...    =   b       .    x   d      ...    +    c      ...  (2-28) 

n-i+j  n-i  n-i+j  n-i+j 


d      ^,'hm-    =    a      .    x   c      ...    +    d      /•,-,>,•,  (2-29) 

n-(i+l)+3  n-i  n-1+3  n-(i+l)+j' 

where  j e[0 , 1, 2 , . . . , i]  are  substituted  in  ascending  order,  and 
(2-28)  is  solved  before  (2-29)  for  each  value  of  j.   Now, 
let  i=2  in  equations  (2-28)  and  (2-29)  and  repeat  the  same 
procedure.   The  index  "i"  is  incremented  until  i=n-l,  and 
(2-28)  and  (2-29)  are  solved  as  before  over  the  appropriate 
range  of  the  index  "j".   The  final  vectors,  C  and  D,  contain 
elements  which  are  the  coefficients  of  the  numerator  and 
denominator  polynomials,  respectively,  of  the  transfer 
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function  (or  driving  point  impedance  function) 


n 

Z      Cs11"1 

•    ,       1 

T(s)     =    =~ (2-30) 

Z      d.sn"j 
j=0      3 


Example 


T(b)  =  10s2+171s+360  (2_31) 

s  +71s  +702S+720 


By  continued  fraction  division: 


10s2  +  171s+360  |s3  +  71s2  +  702s+720^  .1: 


s3+17.1s2+36s 

53.8s2+666s+720 


53.9s2  +  666s  +  720   |l0s2  +  171s+360^   ~~   :  .1855 

10s2+123.562s+133.58 
U7.438S+226.42 


47.U38S  +  226.42   fe  3 . 9s2+66 6s  +  72(A   1.1362: 

53.9s2+257.258s 
408.742S+720 


408.742s  +  720       Jl+7  .  U38S  +  2  26  .  k'Z{       .115 

47.438S+82.79 

143.63  (2-32) 
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/408.742s  +  720^ 
'408. 742s 


143.63        /408.742s+720i.       2.8458s 


720        p.4  3.6  3^        .1 
'143.63 


720 


995 


0 

Therefore,  the  transfer  function  in  the  form  of  equation 
(2-3)  is: 

T(s)  =  1 

.Is  +  1 

.1855  +  1 


1.1362s  + 


.115  + 


2.8458s  +   1 


.1995 

(2-32) 


with 


h,  =  .1  h4  =  .115 

h2  =  .1855  h5  =  2.8458 

h0  =  1.1362  hc  =  .1995.  (2-33) 

3  6 


Now,  using  equations  (2-21)  and  (2-22); 
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a.,  =  h,  =  .1  b,  =  h2  =  .1855 

a2  =  h3  =  1.1362         t>2  =  h^  =  .115  (2-34) 

a3  =  h5  =  2.8458  b3  =  hg  =  .1995. 

From  equations  (2-25)  through  (2-27), 


c3  =  b3  =  .1995 

d2  =  a3xc3  =  2.8458( .1995)  =  .568  (2-35) 

d3  =  1' 

Substituting  i=l  into  equations  (2-28)  and  (2-29)  for  j=0: 


c2  =  b2xd2+c2  =  .115C.568)  +  0  =  .0653 


d,  =  a2xc2+d   =  1.1362C .0653)  +  0  =  .0742 


for  j  =  1: 


c3  =  b2cd3+c3  =  .115(1)  +  .1995  =  .3145 


d2  =  a2xc3+d2  =  1.1362C .3145)  +  .568  =  .9353 


At  this  point,  j=i,  therefore  increment  index  "i" 
for  i=2,  j=0: 


c1  =  b1xd1+c   =  .1855(.0742)  +  0  =  .013; 


dQ  =  a1xc1+d0  =  .K.0138)  +  0  =  .0013; 
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for  j=l: 


c2  =  b1xd2+C2  =  .1855C.9353)  +  .0653  =  .2388 


d,  =  a1xc2+d1  =  .K.2388)  +  .0742  =  .0981 


for  j=2: 


c3  =  b1xd3+c3  =  .1855(1)  +  .3145  =  .5000 


d2  =  a,xd  +d   =  .K.5)  +  .9353  =  .9853 


Now,  at  the  point  where  j=i,  and  i=n-l,  the  transfer 
function  is: 


TCs)  =    .0138s2+.2388s+.5 (2_36) 

,00138s  +. 0981s  +.9853S+1 


Multiplying  numerator  and  denominator  by  1/d   =  7  20  yields: 

T(s)  a   10s2+171s+360  (2_37) 

s  +71s  +702s+720 


2 .   Inversion  of  Cauer  Second  Form 

Initialize  the  two  (n+lxl)  vectors  C  and  D,  (2-23) 
and  (2-24),  to  all  zeros,  except: 


c   =  b  (2-38) 

n    n 

d  ,  =  1  (2-39) 

n-1 

d   =  a  xc  (2-40) 

n    n  n 
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The  following  set  of  equations  are  first  solved  for  i=l. 


c   ....  =  b   .  x  d   ...  +  c  .....  (2-i+D 

n-l+j    n=i    n-i+j    n-i+j+1 

d      r-.T\j_'    =   a       .    x   c      /.,-,>,•    +    d       ...  (2-42) 

n-d+D+j  n-i  n-d+D+j  n-i  +  1 

where  jE[0 ,1 , 2 ,...., i)   are  substituted  in  ascending  order, 
and  (2-41)  is  solved  before  (2-42)  for  each  value  of  j  . 


Next,  find  d   according  to: 
'       n         & 

d   =  a   .  x  c  .  (2-43) 

n    n-i    n 

Now,  let  i=2  in  equations  (2-41)  and  (2-42)  and  repeat  the 
same  recursive  procedure.   The  index  "i"  is  incremented  by 
one  until  i  =  n-1,  and  for  each  value  of  i,  (2-41)  and  (2-42) 
are  iteratively  solved  over  the  appropriate  range  of  the 
index  "j".   The  resulting  elements  of  C  and  D  are  the 
coefficients  of  the  numerator  and  denominator  polynomials 

of  the  transfer  function: 

n 

.,     n-i 

E   c  s 

i 


T(S)  =  hzl (2-30) 

n 

Z      d.i 
3  =  0   1 


n 

_     n-1 

E   d.s 


Using  the  sample  example,  (2-31); 


„,  v     10s2+171s+360 
TC  s ;  = 


s3+71s2+702s+720 


+  .    .  . 

for  3=1,  c  •.-,-,         n 
J   '   n-i+j+1  =  0. 
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Place  the  numerator  and  denominator  terms  in  ascending 
order,  invert,  and  perform  continued  fraction  division 


360  +  171s  +  10s2   /720+702s+71s2  +  s3l-   2 

720+342s+20s2 


360s+51s2+s3 


360s+51s+s3   /360+171s+10s2^   1/s 

2 


60+51s+s 


120s+9s2 


120s  +  9s2   /360s  +  51s2  +  s3^   3 

360s+27s2 

2   3 
24s  +sJ 

24s+s2   /120+9sC   5/s  (2-44) 

120+5s 
45 


4   /24  +  sC   6 
24 


s   M   *+/s 

4 
0 


The  transfer  function  (2-31),  in  the  form  of  equation  (2-6) 
is  : 

T(s)  =  1 

2  +  s 

1+  s 

3+  s 
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5+  

6  + 


(2-45) 


where 

h1  =  2  h4  =  5 

h2  =  1  h5  =  6 

h3  =  3  h6  =  4  (2-46) 

For  the  inversion  process,  using  equations  (2-21),  (2-22) 
and  (2-46) ; 


a1  =  h1  =  2  t^  =  h2  =  1 

a«  =  h~  =  3  b^  =  h^  =  5 

a0  =  h,  =  6  b  =  h  =  4  (2-47) 

3     b  3  b 


and  from  equations  (2-38)  through  (2-40), 

G3  =  b3  =  4 
d2  =  1 

d3  =  a3  x  C3  =  6(4)  =  24.  (2-48) 

Now,  substituting  i  =  1  into  equation  (2-41)  and  (2-42),  for 
j  =  0: 


c2  =  b2  +  c3  =  5  +  4  =  9 

dl  =  a2  x  cl  +  d2  =  3(0)  +  1    =  1  (2-49) 
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for  j=l; 


c:3  =  b2  x  d3  +  0  =  5(24)  =  120 


d2  =  a2  x  C2  +  d3  =  3(9)  +  24  =  51 


and  from  equation  (2-43): 


d3  =  a2  x  c3  =  3(120)  =  360 


for  i=2,  j=0: 


e,  =  b,  =  c2  =  1  +  9  =  10 


dQ  =  a1  x  cQ  +  d1  =  2(0)  +1=1 


for  j=l; 


c2  =  b,  x  d2  +  c3  =  1(51)  +  120  =  171 


d1  =  a   x  c1  +  d2  =  2(10  +  51  =  71 


for  j=2; 


d3  =  b1   x  d3  +  0  =  1(360)  +  0  =  360 


d2  =  a,  x  c2  +  d3  =  2(171)  +  360  =  702 


and  from  equation  (2-43): 


d3  =  a,  x  c3  =  2(360)  =  720. 


(2-50) 
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Therefore , 

C  =  [0    10    171    360] 

D  =  [1    71    702    720],  (2-51) 

and  the  transfer  function  realized  from  equation  (2-30)  is: 

T(s)  -.   10s2  +  171s  +  360  (2_52) 

s  +71s  +702S+720 

which  is  the  same  as  (2-31). 

This  completes  the  development  of  the  continued  fraction 

inversion  algorithms  from  Cauer  First  and  Second  Forms.   This 

iterative  procedure  is  easily  seen  to  be  computationally 

much  simpler  than  Chen  and  Shieh's  method.   First,  it  does 

not  require  the  need  to  find  the  H  matrix,  (2-11);  and 

second,  it  does  not  necessitate  finding  the  coefficients 

of  n  characteristic  polynomials  of  diminishing  order.   This 

method  is  solely  based  on  equation  (2-6),  enumerating  the 

inversion  from  bottom  to  top.   As  by-product,  the  entire 

Routh  array  appears  in  the  intermediate  steps  as  can  be 

seen  from  the  Cauer  Second  Form  example: 


d3 

d2 

dl 

do 

720 

702 

71 

C3 

C2 

°1 

360 

171 

10 

d3 

d2 

dl 

= 

360 

51 

1 

c3 

c2 

120 

9 

d3 

d2 

24 

1 

°3 
1 

4 
1 

(2-53) 
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where  rows  2(n-i)l  and  2(n-i)  are  taken  from  the  ith  iteration, 
i£[0 ,1 , 2 , . . . ,n-l] .   "i=0"  implies  the  rows  come  from  the 
initialization  of  C  and  D.   The  last  row,  the  (2n+l)th,  is 
always  the  single  element  one. 

If  a  comparison  is  made  between  equations  (2-20)  and 
(2-30),  it  is  observed  that: 

6i  =  r2,n-i+l'  i£tl,2,...,n] 

d.  =  rljX_j+1>  je[0,l,2,. .. ,n]  ,  (2-54) 

where  the  c.'s  and  d.'s  are  taken  from  the  (n-l)th  intera- 
tion  under  the  index  "i". 

It  is  also  observed  that  if  the  quotients,  h.'s, 
resulting  from  expansion  into  Cauer  Second  (First)  Form 
are  used  in  the  inversion  algorithm  presented  for  Cauer 
First  (Second)  Form,  then  the  c.'s  and  d.'s  in  the  (n-l)th 
iteration  represent  the  transfer  function  coefficients  in 
reverse  order.   This  is  shown  using  the  preceding  example. 

From  equation  (2-46); 

\   =  2  h4  =  5 

h2  =  1  h5  =  6 

h3  =  3  h6  =  4 

and  equation  (2-47); 

a,  =  h,  =  2      b,  =  h«  =  1 

a2  =  h3  =  3      b2  =  h^  =  5 
a3  =  h5  =  6       b3  =  h6  =  4   . 

31 


Now,  using  the  inversion  scheme  for  Cauer  First  Form,  from 
(2-25)  through  (2-27); 


S3  =  b3  =  4 

d2  =  a3xc3  =  6(4)  =  24  (2-55) 

d3  =  1. 

Making  the  substitution,  i=l,  in  equations  (2-28)  and  (2-29), 
for  j  =  0 ; 


c2  =  t>2  x  d2  +  C2  =  5(24)  +  0  =  120 


d   =  a2  x  C2  +  d,  =  3(120)  +  0  =  360 


for  j=l,  (j=i) 


c'3  =  b2  c  d3  +  c:3  =  5(1)  +4  =  9 


d2  =  a2  x  c3  +  d2  =  3(9)  +  24  =  51 


for  i=2,  j=0; 


c1  =  b,  x  d,  +  c1  =  1(360)  +  0  =  360 


d0  =  al  x  Cl   +   d0  =  2(360)  +  °  =  720 


for  j=l; 


c2  =  b1  x  d2  +  d2  =  1(51)  +  120  =  171 


dx  =  ax  x  c2  +  d,  =  2(171  +  360  =  702 
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for  j  =  2,  (j=i); 

c3  =  b,  x  d3  +  c3  =  1(1)  +  9  =  10 

d2  =  a,  x  c3  +  d2  =  2(10)  +  51  =  71 

and  d~  remains  unchanged,  equal  to  1. 
Therefore; 

•C  =  [0       360       171       10] 

D  =  [720     702        71        1] ,  (2-56) 

and  the  transfer  function  should  be; 


„,,  v      360s2+171s+10 
TC  s)  =  — 


720s3+702s2+71s+l  .  (2-57) 


Since  the  h.'s  from  the  Cauer  Second  Form  were  used  in  the 

l 

inversion  algorithm  from  the  Cauer  Frist  Form,  the  vectors 
C  and  D  require  reversing  non-zero  elements,  resulting  in: 

C  =  [0     10     171     360] 

D  =  [1     71     702     720]  ,  (2-51) 

and  the  correct  transfer  function  is 


T(s)  =    10s2^171s+360  (2.31) 

s  +71s  +7Q2S+720 
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A  digital  computer  program  (FORTRAN  IV)  has  been 
written  for  both  Cauer  First  and  Cauer  Second  Forms  and  is 
included  as  Appendix  3  with  documentation. 
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III.   LINEAR  SYSTEM  ORDER  REDUCTION  VIA  PARTITION 
OF  THE  CAUER  SECOND  FORM 

A  control  system, in  general,  can  consist  of  many  tens 
or  hundreds  of  elements.   In  such  cases,  the  problems  facing 
the  engineer  include:   (1)  too  many  variables  to  efficiently 
handle;  (2)  the  dimension  of  the  system  is  too  high  to  com- 
prehend; and  (3)  the  modifications  needed  to  meet  required 
design  characteristics  are  difficult  to  ascertain.   A  logi- 
cal approach  is  to  seek  procedural  methods  which  reduce  the 
order  of  the  system  to  a  manageable  size  yet  maintain  the 
basic  characteristics  of  the  full  dimension  model. 

A  number  of  different  methods  for  system  simplification 
have  been  proposed  for  the  reduction  of  high  order  dynamic 
systems  to  low  order  models  of  a  more  computationally  or 
analytically  tractable  nature.   The  approaches  used  are 
quite  different,  but  appear  to  fall  into  three  main  groups. 
The  first  is  to  ignore  those  modes  of  the  original  system 
which  contribute  little  to  the  overall  response.   Davison 
[6]  chose  to  neglect  eigenvalues  of  the  original  system 
which  are  farthest  from  the  origin,  retaining  only  the  do- 
minant eigenvalues  and  hence  dominant  time  constants  in  the 
reduced  model.   The  shortcoming  of  this  technique  is  that 
many  systems  do  not  have  any  "dominant"  roots  [7]. 
Chidambara  [8]  essentially  finds  a  reduced  forcing  function 
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so  that  the  steady  state  values  of  the  lower  order  model 
agree  with  those  of  the  original  system.   The  consequence  of 
this  method  is  that  the  approximate  model  give  correct 
steady  state  values  but  incorrect  time  responses  because 
the  reduced  forcing  function  does  not  excite  the  modes  of 
the  two  systems  in  the  same  proportions  [6].   Marshall  [9] 
proposed  the  reduction  of  the  state  matrix  by  partitioning 
it  and  setting  certain  rate  variables  equal  to  zero  in  order 
to  maintain  the  original  steady  state  values.   This  techni- 
que, like  Davison's,  is  based  on  dominant  roots  and,  there- 
fore, exhibits  the  same  shortcomings. 

The  second  approach  is  to  search  in  some  manner  for  the 
coefficients  of  a  set  of  differential  equations  of  specified 
order,  the  response  of  which  is  sufficiently  close  to  that 
of  the  original  system  when  both  are  driven  by  the  same 
inputs.   Sinha  and  Pille  [10]  proposed  a  reduction  technique 
based  on  the  iterative  application  of  the  matrix  pseudo 
inverse  algorithm  [11]  to  determine  a  model  of  specified 
order  which  minimizes  the  sum  of  the  squares  of  the  errors 
between  the  responses  of  the  original  system  and  the  reduced 
order  model  to  a  given  input.   The  main  drawback  of  this 
method  is  that  the  objective  function  to  be  minimized  is 
restricted  to  be  the  sum  of  the  squares  of  the  errors .   Sinha 
and  Bereznai  [12]  presented  a  method  which  minimizes  a 
specified  error  criterion  for  a  given  reduced  order  model 
of  the  original  system,  based  on  the  pattern- search  algorithm 
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of  Hooke  and  Jeeves  [13].   Although  this  method  provides 
more  flexibility  than  that  of  Sinha  and  Pille,  it  generally 
requires  considerably  more  computational  time  due  to  the 
poor  convergence  properties  of  the  pattern-search  algorithm. 

The  third  category  involves  application  of  the  theory 
of  continued  fractions.   Methods  involving  this  approach 
have  been  developed  by  Chen  and  Shieh  [14]. 

Sinha  and  DeBruin  [15]  and  Fellows  et  al  [16]  have 
established  the  fact  that  among  the  methods  previously 
mentioned,  the  approach  by  continued  fraction  expansion  is 
generally  the  best  for  linear  model  simplification. 

A.   SIMPLIFYING  A  TRANSFER  FUNCTION 

The  general  nature  of  a  control  system  is  that  of  a  low 
pass  filter.   Therefore,  model  simplification  should  con- 
centrate on  the  steady  state  aspects  of  the  response  with 
the  transient  portion  given  secondary  consideration.   As 
previously  shown  in  Chapter  II,  the  Cauer  Second  Form 
exactly  characterizes  these  miens. 

Given  the  nth  order  original  system  transfer  function: 

n-1     . 

E    b.s 

T(s)  =  i  =  0   n^ —   ,  (3-1) 

sn  +   E    a.s11 
j  =  0    ^ 

where  an  mth  order  simplified  model  of  the  system  (where 
m  is  strictly  less  than  n)  is  desired,  the  polynomials  in 
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equation  (3-1)  are  rewritten  in  ascending  order; 

■u  j.-u  -L     n-2 .,     n-1 

bn+b1s+....    ....b   0s    +b   ns 

T(s)  =   0   1 x-2 n-1 

an+ans+ a   0sx"2+a   nsn_1+sn      (3"2) 

0—  1  x-z      n-l 

and  expanded  into  a  continued  fraction: 
T(s)  =  1 


hl  + 


h2  ♦ 


h3  + 


h2n-i  + 


h2n  (3-3) 


An  mth  order  simplified  model  is  obtained  by  keeping  the 
first  2m  quotients  of  the  expansion,  omitting  the  remainder; 

T(s)  =  1 


hl  + 


h2  + 


h2m-l  + 


h2m  (3-U) 


and  performing  the  inversion  of  the  truncated  fraction.   The 

inversion  technique  presented  in  Chapter  II  can  be  used  for 

this  purpose. 

3  8 


Consider  the  seventh-order  system,  representing  the 
control  system  of  the  pitch  rate  of  a  supersonic  transport 
aircraft  [10],  described  by  its  transfer  function: 


„,    s  375000Cs+. 08333) 

T(  s;  - 


s7+8  3.6  3  5s6+40  9  7s5+70  34  2s4+8  5  3  70  3s3 


+  2814271s2  +  3310875s+281250  (3-5) 


By  continued  fraction  expansion: 
T(s)  =  1 


9.00036  + 


-.486286  + 


-.036856  + 


781496.2032  + 


.00071478 


(3-6) 
Suppose  a  second-order  simplified  model  is  desired.   Equa- 
tion (3-6)  has  fourteen  quotients.   For  the  desired  system, 
the  first  four  quotients  are  kept  with  all  others  discarded. 
The  truncated  continued  fraction  becomes : 

T^Cs)  =  1 

9.00036  +  s 


-.486286  +      s  (3-7) 


-.036856  +     s 

.616185 
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and  converted  into  transfer  function  form; 


T1(s)  =  — 


.1299s+. 01105 


s+1.14644s+. 09941 


(3-8) 


The  block  diagrams  of  equations  (3-5)  and  (3-8)  in  the 
Cauer  Second  Form  are  shown  in  Figures  3.1  and  3.2,  respec- 
tively.  The  unit  step  and  impulse  responses  of  the  original 
and  simplified  systems  are  compared  and  shown  in  Figures 
3 . 3  and  3.4. 


B.   STATE  EQUATION  SIMPLIFICATION 

The  method  of  system  simplification  just  presented  is 
especially  advantageous  when  converted  into  state  space 
form.   In  Figure  2 . 3,  a  name  for  each  state  variable  is 
given  after  each  integrator,  shown  in  Figure  3.5,  from 
which  the  state  equations  and  output  equation  can  be 
directly  written. 
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C  =  Ch0  h.,  h_....h0  ] 
I      4   b      2n 


'1 


n 


(3-10) 


z  =  Hz  +  Dr  ,  and  C  =  Lz  . 


(3-11) 


Simplification  of  the  equations  in  (3-11)  can  be  achieved 
by  partitioning  of  H,  D  and  L,  as  indicated  in  Figure  3.6. 
The  resulting  mth  order  system  becomes: 


Z"    =   Hp   z      +    D  r    , 
~p        ~F   ~p        ~p 


(3-12) 


where : 


Hp   = 


h2hl 


h4hl 


h2mhl 


h0h,         h.,(h,+hQ)    h0m(h1+hQ) 

21  413  2m      13 


h2hl        VW     h2m(V"--^2m-l) 


(3-13) 
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D" 


(3-14) 


and  C«=  h  z~ , 

p    ~p  ~p» 


c  i  =  [h0  hk 
p       I       4 


h2m] 


pl 

P2 


pm 


(3-15) 


(3-16) 


As  an  example,  consider  the  seventh  order  system  described 
by  the  transfer  function: 


C(s)  .  14i4l.  53s  3  +  783 19  s2  +  525286.125s  +  607693.  25 
R(s)      s7+112.04s6  +  3755.92s5+39736.73s'4 


(3-17) 


+  36  36  50.56s3+7  59894.19s2+68  36  5  6.2  5s+61749  7.3  75 


Arranging  the  polynomials  in  ascending  order  and  expanding 
into  the  continued  fraction  yields : 
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C(s) 
RTsT  " 


1.016133  + 


4.054112  + 


-.067134  + 


(3-18) 


595660.646  + 


.0000757 


5 


and  equations  (3-9)  and  (3-10)  are  formulated  from  the 
quotients.   A  simplified  model  of  second  order  is  desired, 
therefore,  the  state  and  output  equations  are  partitioned 
as  indicated  in  Figure  3.7.   The  simplified  transfer 
function  is: 


Cp(s)     .250367s+l. 035264  , _  no. 

ppf-S  -  —9 (.3-19; 

*rvb;    s+.509768s+l. 051966 


Unit  step  and  impulse  responses  of  the  original  seventh  and 
simplified  second  order  systems  are  shown  in  Figures  3.8  and 
3.9  respectively.   It  is  observed  that  the  results  of  ex- 
pressing the  seventh-order  system  by  a  second  order  model 
are  satisfactory. 

It  should  be  pointed  out  that  a  stable  transfer  func- 
tion may  produce  an  unstable  simplified  function  because  the 
method  of  continued  fraction  expansion  approximation  does 
not  necessarily  guarantee  a  stable  system. 
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IV.   DESIGN  OF  OPTIMAL  LINEAR  CONTROL  SYSTEMS 
WITH  PRESCRIBED  EIGENVALUES 

A.   INTRODUCTION 

Consider  the  control  of  a  plant  with  dynamics  described 
by  a  set  of  first  order,  time-invariant  linear  differential 
equations  of  the  form 

x  =  Ax  +  Bu,  (4-1) 

where  x  is  the  nth-order  state  vector,  A  is  the  (nxn) 
plant  matrix,  u  is  the  scaler  control  and  B  is  the  (nxl) 
input  matrix.   The  output  is  defined  as 

y  =  Cx  ,  (4-2) 

where  C  is  the  (lxn)  output  matrix. 

A  linear  feedback  control  law  is  assumed,  and  of  the  form 
u  =  G*x.+  (4-3) 


There  are  mainly  two  separate  approaches  in  the  deter- 
mination of  the  feedback  control  matrix,  G*,  corresponding 
to  the  system  under  consideration;  1  -  optimal  control  and 
2  -  modal. 

In  the  optimal  control  approach  a  performance  index  is 
considered  which  is  to  be  minimized  in  the  design  of  a 

All  states  are  available  or  an  observer  or  Kalman  filter 
is  used  to  obtain  the  unknown  states. 
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system.   Assuming  a  performance  index  can  be  defined  that 
represents  most  of  the  design  requirements,  "the  solution 
to  the  optimal  control  problem  can  usually  be  obtained  only 
by  numerical  methods  that  yield  solutions  to  only  a  parti- 
cular problem"  [17].   If  solutions  are  sought  to  more  than 
one  numerical  problem,  simple  performance  indices  must  be 
defined,  which  often  do  not  satisfy  many  of  the  design 
requirements.   Therefore,  the  choice  of  a  performance  index 
must  fall  somewhere  between  a  realistic  criterion  and  one 
that  is  mathematically  tractable. 

A  quadratic  performance  index  will  be  considered  as  a 
criterion  for  designing  linear  systems,  of  the  form 

oo 

J    -    h      f         [xT  Q  x  +  Ru  ]dt  ,  (4-4) 

o 

where  Q  is  a  diagonal  non-negative  definite  (nxn)  matrix, 

and  R  is  a  positive  scalar. 

In  the  modal  approach,  G*  is  chosen  so  that  the  closed 

loop  system  achieves  the  prescribed  eigenvalues.   Equations 

(1)  and  (3)  together  yield 

x  =  (A  +  BG*)x  . 

If  Q  and  R  are  given  in  the  optimal  control  approach, 
then  the  eigenvalues  of  the  closed  loop  system  are  uniquely 
determined,  which  may  not  realize  the  required  performance 
characteristics  or  desired  degree  of  stability  for  the 
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system.   Using  the  modal  control  approach,  a  feedback  control 
matrix  can  be  found  that  will  give  the  system  the  desired 
eigenvalues.   This  matrix  is  usually  not  unique,  and  it  is 
not  possible  in  a  practical  sense  to  find  one  that  is 
"better"  than  it's  predecessors,  since  a  performance  measure 
is  generally  not  known  that  corresponds  to  a  given  feedback 
control  matrix.   Therefore,  it  is  necessary  to  find  a  method 
for  determining  the  matrix,  6*,  that  simultaneously  satisfies 
the  desired  eigenvalues  and  minimizes  a  given  performance 
index. 

In  addressing  this  problem,  Chang  [18]  and  Tyler  and 
Tuteric  [19]  have  applied  the  root  locus  method  to  single- 
input,  single-output  and  multivariable  systems,  respectively. 
The  method  lacks  a  rational  computational  procedure  for 
determining  the  elements  of  the  weighting  matrix,  Q,  to  meet 
a  set  of  prescribed  closed  loop  eigenvalues.   Anderson  and 
Moore  [20]  presented  a  restrictive  method  whereby  a  set  of 
eigenvalues  may  be  located  to  the  left  of  a  line  parallel 
to  the  imaginary  axis  in  the  complex  plane.   Chen  and  Shieh 
[14]  presented  a  method  using  sensitivity  analysis. 
Solheim's  [17]  method  of  a  diagonalized  (decoupled)  system 
becomes  complicated  when  the  system  contains  either  complex 
or  multiple  eigenvalues.   Systems  that  cannot  be  diagonalized 
add  further  to  the  complication  of  the  method. 

The  method  developed  here  takes  advantage  of  the 
properties  of  the  Cauer  Second  Form,  is  approached  in  a 
simplistic  manner,  and  is  easy  to  implement  computationally. 
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B.   TRANSFORMATION  TO  PHASE  VARIABLE  FORM 

Consider  an  nth  order  linear  system  of  the  form 

e  =  Se  +  Tf  , 


(4-5) 


with  output 


d  =  We 


(4-6) 


Silverman  [21],  et  al . ,  have  shown  that  if  the  system  is 
controllable,  then  there  exists  a  non-singular  transformation 
matrix  which  takes  an  arbitrary  state  variable  system  to 
phase  variable  (canonical)  form.   (see  Appendix  A) 


x  =  Ax  +  Bu 


(4-1) 


y  =  C  x    , 


(4-2) 


where 


A  = 


0     1     0 0 

0      0      1      0 0 

0      1 

•  •  •  • 

•  •  •  • 

0 
0      0      0 1 

a  ,    a  0    a  0 oc 

nl     n2     n3  nn 
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B  = 


C    --    m±B2     63 


6n] 


If  the  system  is  not  controllable,  the  phase  variable 
(canonical)  form  may  still  be  obtained  from  the  system 
transfer  function 


Y(s)  _ 

U(s) 


n 
i  =  lSiS 


i-1 


n 

n    _ 
s   +  l      a  .  s 

i=l   nl 


(4-7) 


i-1 


Once  the  system  is  in  phase  variable  form,  Chen  and  Shieh 
[22]  have  shown  that  the  equivalent  system  in  Cauer  Second 
Form  is  easily  written  as 


z  =  Hz  +  Vu 


(4-8) 


y  =  C*z  , 


(4-9) 


where  the  two  forms  are  related  by  a  linear  nonsingular  trans^ 
formation  matrix,  P, 


z  =  Px  . 


The  matrix  P  is  an  (nxn)  upper  triangular  matrix. 


(4-10) 
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The  performance  measure  under  consideration  becomes 

oo 

J  =  %      f         [zTQz  +  uTRu]dt  (4-11) 


where 


-IT     -1 

Q  =  (P    )    ()  P 


From  optimal  control  theory,  the  Hamiltonian 

H  =  ^[zTQz  +  uTRu]  +  PT[Hz  +  Vu]  ,  (4-12) 


where  P  is  the  set  of  Lagrange  Multipliers  (also  called  the 
costate  or  adjoint  vector) .   For  the  Hamiltonian  to  be 
globally  minimized,  assuming  no  bounds  on  admissible  states 
and  control  values,  it  is  necessary  that  dH/d u  =  0  and 
32H/au2>0. 

3H/3u  =  Ru  +  VTp  =  0                               (4-13) 

implies 

u*  =  -R_1VTp  ,  and  (4-14) 

92H/3u2  =  R>0  ,  (4-15) 

since  R  was  defined  as  a  positive  scalar.   Included  in  the 
necessary  conditions  for  optimality  are 

z    =    Hz    +    Vu*  (4-16) 

8H/9z    =    -p    =    Qz    +    HTp       .  (4-17) 
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Combining  equations  (4-14),  (4-16)  and  (4-17)  yields  a  set  of 
2n  linear  differential  equations  forming  the  canonical  sys- 
tem in  Cauer  Second  Form. 


-1  T 
H     VR  ^V1 


-Q     -H 


(4-18) 


It  remains  to  be  shown  that  this  form  is  useful  in  obtaining 
optimal  closed  loop  solutions  that  correspond  to  a  set  of 
prescribed  eigenvalues. 

C.   SIMILAR  EIGENVALUES 

Consider  the  2nth  order  cononical  system  in  phase 
variable  form: 


x 


-Q 


-1  T 

-BR  XB 


-A" 


x 


(4-19) 


It  is  known  that  this  system  possesses  n  eigenvalues  with 
negative  real  parts  and  n  eignevalues  with  positive  real 
parte,  and  that  they  are  located  symmetrically  about  the 
imaginary  axis  in  the  complex  plane  [  17 .']  .   The  eigenvalues 
of  the  optimal  feedback  system 


x  =  (A  +  BG*)x 


(4-20) 
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are  identical  to  those  with  negative  real  parts  in  the 
canonical  system.   It  is  possible,  therefore,  to  focus 
on  the  canonical  system  where  the  dependence  eigenvalues 
on  the  weighting  matrices  Q  and  R  can  be  studied  without 
solving  the  matrix  Riccati  equation.   The  problem  is  in 
determining  a  weighting  matrix,  Q,  such  that  the  system 
attains  the  prescribed  set  of  eigenvalues. 

The  canonical  system  in  Cauer  Second  Form  can  be  ob- 
tained from  the  phase  variable  form  using  the  linear  trans- 
formation 


z  =  P  x 


(4-10) 


where  P  is  a  nonsingular  matrix.   We  have 


-Q 


-IT-.     r-   -i 

-BR  ^B        x 


-A' 


(4-19) 


LP  J 


0 


-1  T 
(P  V 


A 


-Q 


-  1  T 

•BR  "B 


-A 


x 


(4-20) 


,-1 


0 


(4-21) 
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A 

P 


PAP 


-1 


-ITT 
-PBR      B    P 


-CP"1)1,   Q    PT         -(P"1)TATPT 


(4-22) 


H 


-Q 


-1   T 
-VR    ^V 


•H" 


(4-23) 


Let 


M   = 


-Q 


-1   T 
-BR    XB 


-AJ 


-1    T 
H  -VR   XV 


N    = 


-Q 


-H" 


and 


P    = 


0 


(P-X)T 


(4-24) 
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where  each  sub-matrix  of  M,  N,  and  P  are  known  to  be  (nxn) 
square  matrices.   It  is  easily  seen  that  P  is  non-singular, 
and  that 


P  P 


-1 


0 


(p-1)1 


«-l 


0 


0 


P  P"1    0 


0 


Cp-l)TpT 


0 


(4-25) 


Therefore , 


P  M  P"1  =  N 


(4-26) 


shows  the  similarity  of  the  M  and  N  matrices.   Two  similar 
matrices  have  the  following  properties: 
1.   Their  determinants  are  the  same. 


det  M  =  det  N 


(4-27) 
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2.   Their  traces  are  the  same 


Tr  [M]  =  TrCN] 


(4-28) 


k 

E 
i  =  l 


m.  .  =   E   n .  . 
j  =1   JJ 


(4-29) 


3 .   Their  characteristic  equations  are  the  same 


det[AI-M]  =  detCXI-N]  =  0   , 


(4-30) 


where  X    is  an  arbitrary  variable.   Since  their  characteristic 
equations  are  the  same,  the  eigenvalues  of  M  and  N  must  be 
identical.   It  is  now  known  that  the  Cauer  Second  Form  and 
phase  variable  system  matrices  are  similar  in  that  they 
possess  identical  eigenvalues. 


D.   DEVELOPMENT  OF  THE  PRESCRIBED  EIGENVALUE  PROBLEM 

Initially  given  is  a  known  linear  system  with  dynamics 
described  by  either  a  set  of  first  order  differential  equa- 
tions or  its  transfer  function.   It  is  desired  to  find  the 
optimal  feedback  control,  u* ,  such  that  the  performance 
measure 


J  = 


hi  f     C  x 


*11 

0 


0    0 

^22  ° 

0    q 


nn 


x   + 
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u[l]u]  }  dt 


(4-31) 


is  minimized,  where  the  eigenvalues  of  the  optimal  system 
are  specified  as 


A-|3  A«,  A»j   •  ■  •  ,  A 


1.  Evaluation  of  the  State  (H)  and  Linear 
Transformation  (P)  Matrices 

-''■■■"----  '  '-     MJ-    ■  ■ 

Assume  the  transfer  function  of  the  known  system  is 


given 


T(s)  = 


n 

S  &  .s 
i  =  l   x 


i-1 


s   +   E    •  s 


i=l  x 


(4-32) 


The  system  matrices  in  phase  variable  form  are 


A  = 


0 
1 


"al    "a2    "a3 


0 
0 

-a 


n 


B  = 


(4-33) 


C=   L3  -I  8  o  •  •  •  •  S  _  J 


(4-33) 


From  equation  (4-32),  the  Routh  array  is  formed 
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al     a2 


a     1 
n 


Sl     62 


.Sn    o 


c    c 


dl     d2 


el     e2 


(4-34) 


In  matrix  notation,  the  Routh  array  becomes  [r..],  where 
i'eCl,  2,  ,  2(n+l)] 


j:e[i,  2,  ,  n+i]  , 

and  elements 

r2(n-k)+3,  k  =  1 

r2(n-k)-4,  k  =  ° 
ke[l,  2,  .  .  .  . ,  n+1]  . 

In  general,  the  (2L+l)th  row  of  the  Routh  array  is  the  Lth 
row  of  P . 
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p  = 


T         r         T 

r31  32     r33 


r51     r52 
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3  ,n 


5,n-l 


7,n-2 


2n+l,2 

'2n+l,l 


(4-35) 


P  being  an  (nxn)  upper  triangular  matrix,  will  always  have 
an  inverse,  P   ,  which  can  be  quickly  and  efficiently  deter- 
mined by  digital  computer.   The  H  matrix  formulation  becomes 


H  =  P  A  P 


-1 


(4-36) 


hll    h12    h13 


hll    h22    h23 


h22    h33 


hll    h22    h33 


In 


2n 


'3n 


n-l,n 


n,n 


(4-37) 


The  elements  of  the  H  matrix  can  also  be  obtained  more 
easily  and  directly  from  the  first  column  of  the  Routh  array. 
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Let 


1+1,1 


(4-38) 


The  h.'s  correspond  to  the  quotients  of  the  continued  frac- 
tion expansion  of  the  transfer  function  in  (32). 


T(s)  = 


hl  + 


h2  + 


h3  + 


\  + 


(4-39) 


The  H  matrix  then  becomes  as  shown  in  Figure  4.1.   The 
regular  pattern  of  the  elements  enable  the  H  matrix  to  be 
obtained  by  inspection  once  the  h.'s  have  been  determined 
from  either  (4-34)  and  (4-38),  or  (4-39). 

The  matrices  V,  and  C* ,  are  easily  obtained: 


V  =  P  B  = 


(4-40) 
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C*  =  c  P"1  =  [h0   h.   hR  h0  ] 

~  ~        2    4b        2n 


(4-41) 


B.  Evaluation  of  Q 


-IT     -1 

Q  =  (P   )   Q  P 


vll 
Ql3 


:12 


22 


23 


Qnl    Qn2 


13 


23 


33 


■n3 


In 


2n 
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•nn 


(4-42) 


The  canonical  system  in  Cauer  Second  Form  (4-23)  has  now 
been  obtained,  with  the  numerical  values  of  the  elements 
of  Q  still  to  be  found.   This  system  will  be  compared  to 
a  non-optimized  system  with  the  prescribed  eigenvalues  also 
in  Cauer  Second  Form.   The  desired  system  in  phase  variable 
form  is : 


x  =  A*x  +  B*u 


y  =  E  x 


(4-43) 
(4-44) 


Formulation  into  Cauer  Second  Form  yields 


x  =  H*z  +  V*u 


y  =  E*z 


(4-45) 
(4-46) 
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In  a  "nonoptimized"  system  (Q  and  R  set  equal  to  0) ,  the 
canonical  system  appears  as: 


H*     0 


-(H*)T 


(M--47) 


where 


H' 


0 


-(H*)T 


N: 


(4-48) 


possesses  the  n  prescribed  eigenvalues  with  negative  real 
parts  and  n  eigenvalues  with  positive  real  parts,  symmetric 
about  the  imaginary  axis  in  the  complex  plane. 

By  equating  the  characteristic  polynomials  of  N  and  N* , 
(4-23)  and  (4-47)  respectively,  it  is  now  possible  to 
determine  the  elements  of  the  Q  weighting  matrix. 


detCsI-N]  =  detCsI-N*] 


det 


-1  T 
sI-H      VR  -V 


sI+H" 


(4-49) 


det 


sI-H* 


0 


0 


sI+(H*)T 


(4-50) 
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E.   DETERMINATION  OF  THE  OPTIMAL  CONTROL  LAW 

Finding  the  elements  of  the  weighting  matrix,  Q,  is 
obviously  a  non-trivial  matter  for  all  but  the  lowest  order 
systems.   The  method  developed  for  obtaining  these  values 
is  based  upon  a  succession  of  matrix  building  blocks,  which 
are  individually  computationally  simple. 
Starting  with  the  matrix,  H, 


H  =  [hi.]  (4-51) 


define  a  new  matrix,  T, 


T  =  Cti,]  (4-52) 


where 


t.  .  =  h. .  -  h.  .  i  i   j  (4-53) 

13    i]    11  J 

and 

n 
t. •  =    I  t. .    .     i  i   n  (4-54) 


li 


j=i+l    ^ 


The  matrix,  T,  therefore,  is  an  (nxn)  upper  right  triangular 
matrix,  where  the  diagonal  elements  are  equal  to  the  sum 
of  all  other  elements  in  the  same  row.   The  next  "building 
block"  matrix,  G,  is  defined  by: 


G  =  [gi.J  (4-55) 

g.  -.  =  1.0  (4-56) 
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for  ie[l,  2 ,  . . . ,  n]   , 


Si, 2  "  '1-2,1-1  ("-57) 

for  j  til,    2 ,  . . . ,  n-1]  , 


n-j  +  2 
Sin  =    Z     (t-|  v  x  §k  i  1}  (4-58) 


for  je[3,  4,...,  x]  and  iE[l,  2,...,  n-j+1],  where  the 
index  j  is  held  fixed  for  each  summation  over  the  index  i. 
The  matrix,  G,  is  an  (nxn)  upper  left  triangular  matrix 
characterized  by  the  first  column  being  all  ones.   One  more 
matrix  needs  to  be  defined  at  this  point.   Let  the  matrix, 
W,  be  defined  as: 


W  =  [w...  ]  (4-59) 

-  XJK 


where 


"i.i.k  ■  (-1)3  x  «i,j  u-60) 

for  i,  j  and  ke[l,  2,  ...,  n]  . 

The  matrix,  T$,  is  therefore  a  tridimensional  array  with  each 
"level"  an  upper  left  triangular  matrix.   Examples  of  the 
matrices,  T,  G  and  W  are  shown  below.   Although  not  evident 
at  this  point,  the  G  and  W  matrices  will  be  used  heavily  in 
obtaining  the  values  of  the  elements  of  the  Q  matrix. 
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From  the  linear  transformation 


-IT     T 
Q  =  (P  V  Q  P1 


(4-61) 


it  is  observed  that  once  the  element  Q..  is  known,  the 
remaining  elements  in  the  same  row  (column)  can  then  be 
obtained  through  a  process  of 
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G  = 
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(4-63) 
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linear  combinations  of  previously  determined  values.   The 
Q  matrix  is  thus  found  in  the  following  manner: 


-  -  Q1X  -  -  Q12  -  -  Q13  -  -Qm Qin- 


Q22  ~  ~  Q23  "  "  Q24  Q2n" 


Q33  "  "  Q34 Q3n 


'44.  ---------  -Q^^ 


where  for  an  nth  order  system  (n>l) 


2n      9    2n      9 

n  (h.)  -  n  (h.r 
Q     1=1 


for  n  =  l,  Q  x  =  (H2H1)2  -  (h^)2 


^nn 


(4-65) 


■11         2n      ? 

n  (h.) 

i=3 


2n   H 

.n   (TT)2   ^H2H1)2-(h2h1)2]  +  (4-66) 

1=3    i 
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3"1   Pk  n-1     "      + 
Q    =  -  I  *>n  l  Q  (4-67) 

13      k=l   ^n-l.n-l    lk 


for  i^j  or  j/n  , 


Qin  =  -  y         Q±j  (4-67) 

3=1     J 


and 

n  n  ^        n       ?      n      0 

V  "  -    *  *    ly  +  (    *  Hii  >    -  <    *    h3j> 

1  =  1  ]=1    J      1  =  1            3=1   JJ 


n   n  n   n 

+  2C   E    Z   H.-CH..-H..)  -   Z    Z   h . . (ft . . -h. . ) ]   (4-69) 

i=i  3=1  '  1X  D:  i:   l=i  3=1  1:  ~31     1: 


The  values  of  the  diagonal  elements  of  Q  (other  than  Q, ,  ) 
generally  involve  varying  linear  combinations  of  already 
determined  values,  more  easily  expressed  in  terms  of  the  G 
and  W  matrices,  rather  than  the  P  and  H  matrices.   To  aid 
in  determining  the  values  along  the  diagonal  the  following 
labels  are  provided  for  G  and  W. 


V  .  =  p 

1,3 
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(4-70) 
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W,    _  .   W. 
l,n-l,i    l,n,i 


W9  n  1   i    0 

2 ,n-l,i 


7 


(4-71) 
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These  two  arrays  will  provide  the  coefficients  of  each   Q..x 

xi-j 

s  "  necessary.   Assume,  for  example,  that  a  5th  order 
system  is  being  considered  in  equation  (4-50).   The  results 
of  the  expansion  of  both  sides  of  the  equation  results  in: 


n     2i  n     2i 

a   +   E  a.szl  =  b   +    Z   b.s  (4-72) 

°    i=l  x        °     L=l   1 


where  a„  =  1,  b   =  1.   By  equating  a.  ,  and  b.  n  for 

n       n        J      n     &   l-l      l-l 

i.e[l,2,  ...,  n] ,  it  is  possible  to  obtain  Q... 

i.e.        a   =  b  (4-73) 

o    o 

a  and  b   are  the  coefficients  of  s  .   From  (4-70),  it  is 
o      o 

observed  that  the  only  element  in  the  s   column  that  is 
non-zero  is  g,,,  which  appears  in  row  1.   The  1  "indicates" 
that  it  is  necessary  to  only  look  in  "level"  1  of  W,  in  the 
column  corresponding  to  s  .   This  yields  only  the  single 
element  W, r,.   Therefore,  the  coefficient  of  Q, , s   will  be 

sis  x  wisi   •  ("-7k) 

2i 

What  remains  to  be  determined  are  the  coefficients  of  s 

not  involving  the  Q..'s.   Because  of  the  symmetry  of  the 
eigenvalues^  of  both  the  system  being  "optimized"  and  the 
system  with  the  prescribed  eigenvalues,  the  remaining 
coefficients  (those  not  involving  a  Q..)  are  easily 
determined  from: 
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detCsI-H]  x  detCsI  HT]  (4-75) 


and 


detCsI-H*]  x  det[sI+(H*)T],  ("4-76) 


which  give: 


and 


n  n      , 

(   Z    a-  sX)(   Z   (-1)  a.  s1)  (4-77) 

i=0  i=0       1 


n  n      , 

(   E    a.*s1)(   S   (-l)Ka.*s1)  (4-78) 

i=0    x      i=0       x 


where  k=i+l  for  n  even  and  k=i  for  n  odd.   The  indicated 
multiplication  in  (4-77)  and  (4-78)  result  in: 

n      k    2i 

Z   (-1)  a.s  x  (4-79) 

i  =  0        1 
and 

n      k     2i 

Z   (-l)Ka.*sZ1  (4-80) 

i  =  0       x 


respectively,  where  the  same  conditions  are  imposed  on  "k" . 
Returning  to  the  5th  order  example,  the  coefficients  of 
s   are  -an  and  -crn*.   To  obtain  Q,  ,  is  now  a  matter  of 
solving  the  equation: 
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Qll  =  (g    •  g    )  (4-82) 


With  Q, ,  known,  it  is  a  simple  matter  to  obtain  the 
remaining  elements  in  the  first  row  (and  column)  using 

(4-67)  and  (4-68) . 

2     2 
To  find  Q?o>  find  all  the  coefficients  of  s  .   s   results 

in  three  separate  ways: 
M,    0     2 

(1)  S   X  s 

(2)  s1  x  s1  (4-83) 

M,    2     0 

(3)  s   x  s 

The  multiplicand  indicates  which  columns  of  G,  (4-70),  is 
of  interest.  Any  non-zero  element  in  G  tells  which  level 
of  W  is  of  interest.   The  multiplier  is  the  indicator  for 

the  column  of  interest  in  the  array,  W.   Therefore,  for 

0     2 
s   x  s   : 

«15  X  (W131  5X1  +  W231  Q12  +  W331  Q13)  (4-84) 

e  1     1 

for  s   x  s  : 


«1,U  X  (W1U1  Qll  +  W241  Q12>  + 

g    x  (W142  Q21  +  W2,2  Q22)  (U-85) 
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2     0 
for  s   x  s 


«1,S  X  (H151  «U»  + 
§2,3  x  (W152  ^21)  + 
83,3  x  (W153  hi*       •  (U-86) 


Obtaining  Q22  is  now  a  matter  of  solving  the  equation  given 
by  the  combination  of  (4-79),  (4-80),  (-84),  (4-85)  and 
(4-86)  : 


(.Z=1  Sj,3  W15j  Qji5  -  °1  =    "al*  •  (4"87) 

for  Q22. 

Once  Q22  has  been  found,  the  remainder  of  the  2nd  row 

(and  column)  can  be  obtained  using  (4-67)  and  (4-68). 

In  a  like  manner,  the  respective  equations  to  be  solved 
for  Q33,  QU4  and  Q55  are: 


15-        2 
( 
j=i      1=1 


E   g.  -   EW.,.Q..)  +  (Zg.  .,   ZW...  Q..) 
=  1   J'5  i  =  l  ll:   31     i=l  3'   i  =  l  l2D   31 


+  (  Z   «1  3   *  Wi3i  Qii}  +  (  Z  §i  2  -1  Wi4i  Qii} 


5 


+  S=i  gj'x  i=i"i5j  Qji>  "a2 


=  -a*       ,  (4-88) 
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3 
( 

i 


l      g.  .  I   W., .  Q. •)  +  (  I  g.  -   Etf...  Q..) 


+  (.En  gj,l   *Wi3j  Qji)  +  a3  =  a3*  -  (4-89) 

j=l   J  '   i=l    J   J 


(.^gj,i  .^"iij  Oji'  -  %  =  *V    ■  (4-90) 


The  underscore  in  (4-88),  (4-89)  and  (4-90)  indicates  the 
term  that  contains  Q33)  Quii*  anc*  Q55'  respectively. 

The  entire  Q  matrix  is  now  known.   By  the  inverse 
transformation  of  (4-42), 

Q  =  PT  Q  P  .  (4-91) 

Q  along  with  A,  B  and  C  of  (4-1)  and  (4-2)  are  used  to 
obtain  the  matrix  Riccati  equation  steady  state  solution. 

0  =  KA  +  ATK  -  KBR~1BTK  +  Q   .  (4-9  2) 

Once  K  has  been  determined,  the  optimal  control  law  is 
given  by 

u*  =  -R_1BKx   .+  (4-9  3) 

Once  u*  has  been  determined,  an  inverse  non-singular  trans- 
formation can  be  performed  to  take  the  phase  variable  form 

+  It  is  now  known  that  G=-R~1BTK  in  (4-20). 
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back  into  state  variable  form.   (See  Appendix  A) 

F.   ILLUSTRATE  EXAMPLES 

1.   Odd  Order  System  (Third) 
Given  the  system 


—                _ 
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1 

u     (4-94) 


y  =  [1    0    0] 


(4-95) 


find  the  optimal  control  law,  u* ,  such  that  the  quadratic 
performance  index  (4-31)  is  minimized,  where  the  eigenvalues 
of  the  system  are  specified  as 


s,  =  -3,  s2  =  -4,  s3  =  -6. 
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Forming  the  Routrh  array  yields 
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0 

from  which  the  third,  fifth  and  seventh  row  are  extracted  to 
form  P . 


P  = 


19  7  1 
0  30/7  1 
0        0       1 


(4-97) 


and 
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-49/570 

7/30 

0 


1/30 
-7/30 

1 


(4-98) 


From  (4-97)  and  (4-98),  H  and  V  are  calculated 

-13/19  637/570       -13/30 

-13/19  -63/19        9/7 

-13/19  -63/19        -3 


H  =  PAP 


-1 


(4-99) 
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V  =  PB  = 


(4-100) 


The  set  of  prescribed  eigenvalues  yields  the  system: 


x. 


-72    -54    -13 


[~xr 

~0  " 

x2 

+ 

0 

Lx3_ 

1 

u 


(4-101) 


and 


y  =  [1    0    0] 


x. 


(4-102) 


From  the  Routh  array 


72 

54 

13 

1 

1 

0 

0 

0 

54 

13 

1 

-13/54 

-1/54 

0 

115/13 

1 

1/115 

0 

1 

0 

(4-103) 
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extracting  rows  3,  5  and  7  yields  the  P*  matrix 


54  13  1 
0  115/13  1 
0       0       1 


,  (P*)"1 


from  which  H* and  V* are  calculated, 


1/54  -169/6210  1/115" 
0  13/115  -13/115 
0      0  1 

(4-104) 


H*=P*A*(P*) 


-1 


-4/3    676/345      -72/115 
-4/3    -2574/621    1980/1495 
-4/3    -2574/621    -11254/1495 


(4-105) 


V*  =  P*B 


&  DA   - 


1 
1 

1 


(4-106) 


Formulation  of  equation  (4-50)  yields 


det 


si  -  H 


-1  T 
VR  "V 


sI+H" 


=  det 


sI-H- 
0 


sI+H" 
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det 


s+13/19   -637/570   13/30 
13/19     s+63/19    -9/7 
63/19 

Q22 
Q. 


13/19 
^12 


13 


23 


s+3 

Q23 

Q 


l 
l 
l 


l 
l 
l 
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s-13/19  -13/19  -13/19 
637/570  s-63/19  -63/19 
-13/30    9/7       s-3 


=det 


s  +  4/3   -676/345     72/115         0 
4/3     s+2574/621   -1980/1495     0 


4/3 

0 
0 
0 


2574/621 

0 

0 

0 


s+11245/1495   0 


s-4/3 


0 
0 
0 
-4/3 


0 
0 
0 
-4/3 


676/345  g-2574/  -2574/ 

621     621 
-72/115  1980/  s-11245/ 
1493     1495 

(4-107) 


It  is  desired  to  determine  the  values  of  Q.   Brute  force 
enumeration  of  the  determinants  and  equating  coefficients 
of  like  powers  of  s  would  eventually  lead  to  the  solution. 
Using  instead  the  method  developed,  from  H(not  H*)  evolve 
the  T,  G,  and  W  matrices. 


T  = 


7:0000  4.4333  2.5667 
0  4.2857  4.2857 
0       0        0 


(4-108) 
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1  7.0000 
1  4.2857 
1      0 


19.0000 

Q 

0 


(4-109) 


and  W  = 


-1  7.0000 
-1  4.2857 
-1     0 


-19.0000 

0 

0 


(4-110) 


Now  add  the  appropriate  labels  as  in  (4-70)  and  (4-71) 


for  G 


s 

s 

s 

1 

""1 

7.0 

19.0 

2 

1 

4.2857 

0     ' 

3 

1 

0 

0 

(4-111) 


for  W: 


2 

s 

1 
s 

0 
s 

9u 

~-l 

7.0 

-19.0  ~ 

Qi2 

-1 

4.2857 

0 

Si3 

_-l 

0 

0 

(4-112) 


Still  necessary  are  the  coefficients  in  (4-79)  and  (4-80) 


n 


n 


n 


<  I  a.sx)(  z  C-DV.s1)  =   I  (-l)ka.s2j 
i=0  1    i=0      x      i=o      1 


(4-113) 


for  n  odd  (n=3),  k=i+l, 


(  Z  a.s1)(  E  (-l)1+1a.s1)  =   Z  (-l)1+1a.s2i  (4-114) 

i=0      i=0        1      i=0        1 


(s3+7s2+19s+13)(s3-7s2+19s-13) 

=  s6  -  lis4  +  179s2  -  169,  (4-115) 


and 


(  Z  a.^si)(  Z  (-l)i+1a.^si)  =   Z  (-1)  i  +  1a .  *s2i      (4-116) 
i=0       i=0        1       i=0        x 


(s3+13s2+54s+72)(s3-13s2+54s-72) 

=  s6  -  61s4  +  1044s2  -  5184  .  (4-117) 


Now,  by  equating  coefficients  of  s  ,  it  is  possible  to 
obtain  Q -.  -.  -   From  G,  the  only  non-zero  element  in  the  "s 
column"  is  19,  which  corresponds  to  row  1,  and  therefore 
level  1  of  W.   From  level  1  of  W,  the  only  non-zero  element 
in  the  "s   column"  is  -19,  with  row  label  Q1 , .   The 
coefficient  of  Q, ,  as  obtained  from  (4-81)  is  19(-19)  = 
-361.   The  solution  for  Q,,  is  obtained  from  equations  (4-81) 
and  (4-82) : 
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-361  Q    -  169  =  -5184 

Qll  =  (_5184  +  169)  /    ~361 

=  13.892   .  (4-118) 

From  (4-65),  (4-67)  and  (4-68),  Q   and  Q    are  successively 
obtained: 


5l2  =  r^f  *U    =  -WT  (13.892)  =  -22.69  (4-119) 


Q13  =  "Qll  "  Q12  =  -13'892  "  (-22.69)  =  8.798.      (4-120) 


The  next  power  of  s  which  results  from  expansion  of  (4-107) 

2         .         .  .         2 
is  s  .   Equating  coefficients  of  s  ,  it  is  now  possible  to 

"2  0   2    11 

obtain  Q97«   s   results  from  the  products  s  xs  ,  s  xs  , 

2   0 
and  s  xs  .   Therefore,  from  the  development  starting  at 

(4-83),  the  equation  to  be  solved  for  Q22  is: 

j=l  J'    1=1     J    J         j=l  J     1=1     J    J 


+  (  I  g^  ,.   I  W.,.  Q..)    +  a1   =  a,*  (4-121) 


-pji  i  =  1  l33   Ji 
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Substituting  known  values,  equation  (4-121)  becomes: 
Sl3(Wll  +  W211  512  +  W311  313>  + 

§12(W121  Qll  +  W221  Ql2}  +  S22(W122  Q21  *  W222  Q22) 
+  Sll(W131  5U)  +  S21CW132  521>  +  S31(W133  331> 


+  o1   =  ax*  (4-122) 


19.0C-H13.892)  ^  K-22.69)  -  1(8.798)] 


+  7.0  [7.0(13.892)  +  4 . 2857 (-22  .  69  )  ] 


+  4.2857  [7. (K-22.69)  +  4.2857  Q22>]  +  lC-19 (13 . 892 ) ] 


+  1  C-19C-22.69)]  +  1[-19(8.798)]  +  179  =  1044 

Q22  =  84.155  (4-123) 

From  (4-6  8) 

Q23  =  "521  "  Q22  =  -Ql2  "  ^22  = 

-  (-22.69)  -  (84.155)  =  -61.465  (4-124) 
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s   is  the  next  power  of  s  obtained  in  the  expansion  of  (4-107) 

.  .  4 

Equating  coefficients  of  s  ,  it  is  now  possible  to  obtain  the 

4 
equation  from  which  Q~-  can  be  found.   s   results  from  multi- 

2   2 
plying  s  xs  .   Therefore,  from  G  and  W: 


gllC"1(Qll)  '  1(Q12}  '  1(Q13)] 


+  g2it-lCQ2i)  -  1(Q22}  '  1(Q23)] 


+  831C-1(Q31)  "  1(Q32}  "  lC(W 


"  °2    =  ~Q2*   '  (4-125) 


where  g,  ,  ,  g21  and  g-,  all  equal  one.   Therefore,  (4-125) 
becomes 


-   E     E    Q. .  -  a0  =  -a  *  (4-126) 

i=l   j  =  l    ^     2      2 


Solving  (4-126)  for  Q33: 

3      3 
Q„  =  -a  *  +  a      -  I  Z      Q..  (4-127) 

66  l  l        i  =  l    j=l   13 

i+j  t    6 

Since 

k=l 
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n  ^ 

1    Qik  =  °  ' 
k=l  lK 


(4-128) 


and  (4-127)  simplifies  to 


Q3k  =  +  a2*  "  a2  -   J   Q3k 


(4-129) 


=  61-11-8.798+61.465  =  102.667 


(4-130) 


Q  is  now  entirely  known. 


Q  = 


13.892 
-22.69 
8.798 


From  (4-91): 


-22.69 
84.155 
-61.465 


8.798 

-61.465 

102.667 


(4-131) 


Q  =  P  QP  = 


5015 

0 

0 

0 

865 

0 

0 

0 

50 

(4-132) 


To  find  the  optimal  feedback  control  law,  u* ,  it  is  necessary 
to  solve  the  matrix  Riccati  equation  (4-92),  where  Q,  A,  B 

r*  *«      *»■» 

and  C  are  as  given  in  (4-131),  (4-94)  and  (4-95).   The 
solution  yields: 
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u*  =  -[59    35    6] 


(4-133) 


x, 


or  u*  =  -59  x,  -  35  x0    -  6  x~ 


With  G*  as  given  in  (4-133),  x  =  (A  +  BG*)x 
becomes 


xl 

. 

x2 

*™ 

. 

_X3_ 

0  10 
0  0  1 
72    -54    -13 


~xl~ 

x2 

_X3_ 

(4-134) 


Comparing  (4-134)  with  (4-101),  it  is  observed  that  (A+BG*)= 
A:':;  therefore,  the  desired  system  has  been  realized  with  the 
set  of  prescribed  eigenvalues. 

2 .   Even  Order  Example  (Fourth) 

Consider  the  linear  system  represented  by  the  transfer 
function: 


Vs)  = 


s  +9s+34 


s4+12s3+48s2+80s+4 


(4.135) 


the  system  eigenvalues  are  -2,-2,-2,  and  -6. 
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This  system  could  be  expressed  in  phase  variable  form  (4-31), 
thereby  obtaining  the  transformation  matrix,  P;  the  system 
can  be  realized  in  Cauer  Second  Form  (4-36).   Instead,  in 
this  example,  it  is  desired  to  obtain  the  quotients  which 
result  from  partial  fraction  expansion  of  the  transfer  func- 
tion.  In  Cauer  Second  Form,  T(s)  becomes: 


T(s)  =  34±_9s  +  s2 


49  +  80s  +  t8s  +12s3  +  s't 


hl  + 


V 


V 


\   +  2. 


h5  + 


h6  + 


h?  ♦ 


h8 
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(t-136) 


where  the  quotients,  h.,  are 

1 


h- 


1.41176 
.505245 
-4.6287 
-.627911 


•8 


23.07 
.146914 
-281.808 
-.024241 


Substituting  these  values  into  the  matrix  in  Figure  4.1 
yields 


H  =  - 


71328 
71328 
71328 
71328 


-.88647 
2.01997 
2.01997 
2.01997 


.20741 
-.47261 
2.91670 
2.91670 


-.03422 
.07798 
-.48125 
6.35004 


(4-137) 


Again,  V  is  (and  in  all  cases  considered  will  be)  a  vector 
of  all  ones , 

"1 

1 
1 

1 


V  = 


(4-138) 


It  is  desired,  as  before,  to  find  the  optimal  control  law, 
u« ,  such  that  the  quadratic  performance  index  (4-31)  is 
minimized,  and  that  the  optimal  system  realizes  a  set  of 
prescribed  eigenvalues.   The  eigenvalues  are  given  as: 
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2,  s2  =  -5 


s3»  s4  =  ~6  1  J3 


(4-139) 


With  the  zeroes  of  the  transfer  function  (4-136)  the  same, 
the  transfer  function  of  the  desired  system  with  the  pre- 
scribed eigenvalues  becomes: 


V(8)  =  — 


s  +9s+34 


s  '+19s3+13  8s2+435s+45  0 


(4-140) 


By  continued  fraction  expansion  of  T*(s)  the  quotients 
obtained  are: 


H. 


H, 


H, 


H, 


13.2353 

.107635 

-71.3206 

-.088175 


H, 


H, 


H. 


H, 


-1077.43 
.004834 
396.926 
-.024294 


from  which  H*  is: 


H*  =  - 


1.42458  -1.16702  .06399  -.32154 

1.42458  5.12168  -.28082  1.41114 

1.42458  5.12168  -5.48975  27.58655 

1.4258  5.12168  -5.48975  17,94349 


(4-141) 
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Formulation  of  equation  (4-5  0)  yields  in  matrix  notation: 


det  = 


sI-H 
Q 


-.  + 


sI+H" 


det 


sI-H* 

0 


sI+(H*)T 


(4-142) 


which  is  to  be  solved  for  Q. 

As  before,  we  formulate  the  T,  G  and  W  matrices 
according  to  equations  (4-53)  and  (4-54);  (4-56),  (4-57) 
and  (4-58);  and  (4-60),  respectively. 

This  yields: 


T  = 


12.00 
0 
0 
0 


2.9064 

9.6614 

0 

0 


2.7093 
3.3893 
6.8313 
0 


6.3500 
6.2721 
6.8313 
0 


(4-143) 


G  = 


1 

12.0000 

46.5882 

67.2941 

1 

9.6614 

23.1534 

0 

1 

6.8313 

0 

0 

1 

0 

0 

0 

(4-144) 


the  notation  1  indicates  an  (nxn)  matrix  with  all  elements 
unity . 
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w  = 


-1  12.0000 

-1  9.6614 

-1  6.8313 

-1  0 


-46.5882 

-23.1534 

0 

0 


67.2941 

0 

0 

0 


(4-145) 


With  the  proper  labels,  (4-70)  and  (4-71), 
for  G: 


1 
2 
3 
4 


s 
1 
1 
1 
1 


12.0000 
9.6613 
6.8313 
0 


46.5882 

23.1534 

0 

0 


0 
s 

67.2941 

0 
0 
0 


(4-146) 


for  W: 


Qi2 

Qi3 


14   l_ 


,t      s 

"-1  12.0000 

-1  9.6613 

-1  6.8313 

-1  0 


46.5882 
23.1534 
0 

0 


67.2941 

0 

0 

0 


(4-147) 


2k 
G  and  W  contain  all  of  the  coefficients  of  Q..s    resulting 

lj 

from  the  expansion  of  (4-142).   Still  to  be  found  are  the 
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2k 
coefficients  of  s    not  involving  any  Q..,  the  a-'s  and  a  •  * '  s  , 

Therefore,  from  (4-79)  and  (4-80): 


(s4+12s3+4  8s2+8  0s+48)(s4-12s3+48s2-80s+4  8) 


=  s8  -  48s6  -  480s4  -  1792s2  +  2304,  (4-148) 


an  =  2304 

a3  =  -48 

a   =  -1792 
1  a^  =  1  (4-149) 

a2  =  480 
and 

(s4+19s3+138s2+435s+450)   x 
(s4-19s3+13  8s2-43  5s+450) 


=  s8-85s6  +  3414s4  -  65025s2  +  202500   ,  (4-150) 


oQ*  =  202500 


ax"  =  -55025 


a2*  =  3414 


a3*  =  -85 


a4*  =  1       .  (4-151) 
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Determination  of  the  Q..'s  results  from  equating 
coefficients  of  like  powers  of  s  from  the  expansion  of  the 
determinants  in  equation  (4-142).   The  Q..'s  are  obtained 
in  a  successive  method  (4-65)  by  starting  with  the  s 

coefficients,  then  continuing  by  equating  coefficients  of 

2k 
s   ,  kieLl,  2,  ...  n-1],  in  scandent  fashion. 

For  s°: 

from  G,  the  only  non-zero  element  is  g..,  corresponding  to 

row  1,  therefore,  level  1  in  W.   The  only  non-zero  element  in 

level  1  of  W  in  the  s   column  is  W,u,.   Equating  coefficients 

yields : 


(glu  x  Wm)  Qu  ♦  aQ  =  a0*    .  (H-152) 

Solving: 

n    -    202500-2304        :  (4-153) 

Qll  "  (67.2941)(67.2941)   "    ^-20803  Q4-153J 

and  from  (4-6  7)  and  (4-68), 

Q12  =  Q21    =    -88.95327 


^13  =  Q31    =    48'14819 
Q    =  Q41    =    -3.40294 


(4-154) 
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2 
for  coefficients  of  s 


1        3       ,.       2       2 

Z      g.    Z    W   ,  Q    +  Z    g,3  Z   W.  .  Q.. 
1=1   :H  i  =  l  10   31   i=l  ^  i=l  ^   ^ 


+  Z   g.    Z  W.,.  Q..   +   a   =  +a*    .  (4-155) 

j=l  •>*  L=l  x  :   J1      -1     -1 


Solving : 

Q22  =  296.94152  (4-156) 

and  from  (4-67)  and  (4-68), 
Q23  =  Q32  =  -263.70162 
Q24  =  Q42  =  55'71337  (4-157) 


4 
for  coefficients  of  s  : 


^gj3  J/iij  Qji  +   *«j2  J/wj  V 


4  2-. 

+      Eg..      2   W.0.    Q..    +    a      =    a   *  (4-158) 

.    .6ii    T     ,    i3]    x]i  2  2 

j=l    J       L=l        J       J 
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Solving: 

Q33  =  351.24159  (4-159) 

and  from  (4-68) , 


Q34  =  Q43  =  ~135-68816  (4-160) 


for  coefficients  of  s  : 


Eg..   E  W.,.  Q .  .  +  a  0  =  a  Q  * 
j=l  31  i=1  il3   31     3     3 


Q  is  now  entirely  known 


(4-161) 


4    4 

Q44  =  "   E    E  Qii  +  Q3*  '  a3  (4-162) 

44      j=l  1=1  1=l     6  * 

i+j?*8 


-Q44  =  -  .2  QHj  +  *3*  "  a3  (4"163) 

3=1   J 


Q    =  120.37773  (4-164) 


103 


Q  = 


44.20803  -88.95327     48.14819 

-88.95327  296,94152  -263.70162 

48.14819  -263.70162  351.24159 

-3.40294  55.71337  -135.68816 


-3.40294 
55.71337 
-135.68816 
120.37773 


From  equation  (4-91): 


(4-165) 


200196 

0 

0 

0 

9= 

PiQP    = 

0 
0 

63233 

0 

0 
2934 

0 

0 

0 

0 

0 

37 

(4-166) 


Substituting  matrices  Q,  R,  A,  B  and  C  into  the  Riccati 
equation,  and  solving,  yields  the  optimal  feedback  control 
law ,  u*  ; 


u*  = 


=  -[402    355    90    7] 


~xl~ 

x2 

X3 

_v 

(4-167) 


where  G*  =  -[402    355    90    7]. 

The  optimal  closed  loop  system,  x  =  (A+BG*)x, 
in  phase  variable  form  is: 


(4-168) 
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~'xl 

x2 

- 

= 

X3 

• 

_v 

-450   -435   -138   -19 


rv 

X2 

X3 

_V 

(4-169) 


which  realizes  the  desired  set  of  eigenvalues. 
3.   Higher  Order  Example  (Seventh) 

The  previous  examples  represent  an  odd  and  an  even 
ordered  system,  illustrating  the  minor  differences  in  com- 
putational procedures .   It  is  observed  that  for  low  order 
systems,  the  calculations  can  be  done,  relatively  easily  by 
hand.   Higher  order  systems  require,  laborious  and  tedious 
computations.   Appendix  C  provides  a  digital  computer  program 
which  yields  the  weighting  matrix,  Q,  with  the  only  required 
input  being  the  transfer  functions  of  the  known  and  desired 
eigenvalue  systems. 

The  following  seventh  order  example  utilizes  the 
results  from  the  program  given  in  Appendix  C. 

Consider  the  linear  system  given  by  its  transfer 
function : 


T?(s)  =  — 


249.435788 


'+9.0s6  +  40.4s5  +  116.8sl++2  3  3.6s3  +  3  2  3.2s2 


(4-170) 


+  288.0s  +  128.0 
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It  is  known  that  the  characteristics  of  the  desired  system 
are  such  that  its  resulting  transfer  function  is  given  by: 


T    *(s)    _  2i+9  .1435788 


s7  +  15.4s6  +  101,64s5  +  3  7  2.6  8S1+  +  819.8  9  6s3 


+1082.2  6272s2+793.659328s+249.43  57888 


(4-171) 


The  diagonal  Q  matrix  was  determined  from  the  computer  pro- 
gram yielding: 

Q1-  =  45834.21273428 

Q22  =  89780.21841734 

Q33  =  55970.39786199 

Q^  =  19170.7157376 

Q55  =  3959.33664 

Q66  =  494.9776  (1+_1?2) 

Q7?  =  33.68 

Q,  R  and  the  state  and  output  matrices  representing  the 
transfer  function  in  equation  (4-170)  in  phase  variable  form 
were  substituted  into  the  matrix  Riccati  equation.   The  opti- 
mal feedback  control  law,  u*,  resulting  from  solution  of  the 
Riccati  equation  is: 
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u*  =  -121.435789  x±    -  505.659319  x? 

-759.06269  x3  -  586.29596  x^  -  255.879974  x5 

-61.2399914  Xg  -  6.39999944  x?   ,  (4-173) 

where  the  matrix  G*  is: 

6*  =  -[121.435789  505.649319  759.06269  586.29596 

255.879974  61.2399914  6.39999944].  (4-174) 

The  optimal  closed  loop  system, 

x  =  (A  +  BG-)  x 

y  =  Cx    , 
expressed  as  a  transfer  function  is: 

Y(s)     249.435788 

U* ( S }    s7+15.39999974s6+101. 6399914s5  + 

372. 67997 4s4+ 819. 89596s 3+ 

1082. 26269s 2+ 793. 659319S+ 249. 435789, 

which  realizes  the  desired  transfer  function  with  the  pre- 
scribed eigenvalues. 

Figures  4.2  through  4.13  show  the  impulse  and  step  res- 
ponses of  the  three  previous  examples,  both  before  and 
after  compensation. 
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V.   CONCLUSIONS  AND  DISCUSSION 

A  method  has  been  shown  for  determination  of  the  state 
weighting  matrix  in  order  to  satisfy  a  prescribed  set  of 
eigenvalues  through  phase  variable  state  feedback.   From  a 
strictly  mathematical  viewpoint,  this  technique  requires  only 
a  knowledge  of  matrix  algebra.   Every  attempt  has  been  made 
to  avoid  the  necessity  of  inverting  a  matrix.   The  intro- 
duction of  Chapter  IV  made  known  the  fact  that  previous  devel- 
opments in  this  area  have  suffered  the  main  drawback  of  res- 
triction.  The  author  believes  the  method  presented  here, 
using  Cauer  Second  Form,  overcomes  many  of  these  restrictions 
It  presents  a  rational  computational  procedure  for  determina- 
tion of  the  weighting  matrix,  Q;  the  system  eigenvalues  are 
only  required  to  be  in  the  left  half  of  the  complex  plane  as 
opposed  to  the  left  of  a  line  parallel  to  the  imaginary  axis ; 
and  the  method  is  no  more  complicated  for  multiple  or  complex 
eigenvalues  than  a  system  with  linearly  independent  eigen- 
values (or  eigenvectors). 

It  should  be  noted  that  some  authors  "define"  the 
eigenvalue(s)  of  a  matrix  to  be  only  the  real  root(s)  of  the 
characteristic  equation.   In  the  development  that  has  pre- 
ceded in  this  thesis,  all  roots  are  considered  eigenvalues 
of  the  associated  matrix. 
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The  algorithm  derived  in  Chapter  IV,  is  designed  as  a 
basis  for  future  research.   In  particular,  if  the  designer 
is  working  with  nth  order  systems  where  n  is  relatively 
large,  it  may  be  advantageous  to  look  at  using  mth  order 
simplified  models  (m<n)  of  each  system  by  a  partitioning 
scheme  similar  to  that  in  Chapter  III.   Again,  it  is  em- 
phasized that  reduced  order  models  do  not  necessarily  yield 
stable  systems.   If  the  simplified  system  retains  the  basic 
characteristics  of  the  original  system,  especially  in  steady 
state,  then  this  would  appear  to  be  a  reasonable  approach. 

A  parallel  approach  could  also  be  inventigated  regarding 
multi-input  multi-output  (MIMO)  systems,  treating  each 
element  of  the  system  transfer  matrix  as  an  individual 
transfer  function. 

To  the  author's  knowledge,  no  work  has  been  done  in  the 
digital  or  sampled  data  areas  involving  continued  fraction 
theory.   This  area  should  be  considered  due  to  the  increasing 
use  and  need  for  digital  control  systems . 

These  topics  represent  just  a  few  of  the  areas  available 
for  future  work. 
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APPENDIX  A 


Consider  the  system 


X  =  Ax  +  Bu, 


(A-l) 


where  x  is  an  n-dimensional  state  vector,  u  is  the  input 
function,  and  A  and  B  are  time-invariant  (nxn)  and  (nxl) 
matrices,  respectively.   The  phase  variable  (canonical) 
system  representation  is  defined  as 


v  =  Av  +  Bu, 


(A-2) 


where  v  is  an  n-dimensional  state  vector  and 


A  = 


0  10 
0  0  1 
0    0    0 


"*1  ~a2  ~a3 


n 


,   B  = 


(A-3) 


The  systems  represented  in  (A-l)  and  (A-2)  are  said  to  be 
equivalent  if  and  only  if  there  exists  a  non-singular  matrix, 
K,  such  that 


x  =  Kv 


(A-4) 
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Kalman  [23]  has  shown  that  a  necessary  and  sufficient  condi- 
tion for  such  an  equivalence  to  exist  is  that  the  system 
in  (A-l)  be  completely  controllable. 

The  controllability  matrix  of  system  (A-l)  is  defined  by 


E  =  [B  :  AB  I  A  B 


A^B] 


(A-5) 


or   in  an   equivalent  manner 


hi    —    L  e  -,       e  ~       .  .  .  .       6    j  j 

12  n 


(A-6) 


where  the  (nxl)  vector  e.  is  recursively  defined  as 


~i+l   ~  ~i 


e,  =  B  . 


(A-7) 


The  controllability  matrix  of  system  (A-2),  E,  is  defined 
in  a  similar  manner  with  A  and  B.   Since  there  is  only  one 
control  input,  a  necessary  and  sufficient  condition  for 
controllability  is  that  the  (nxn)  matrix  E  (or  E)  have  an 
inverse. 

Silverman  [20]  has  shown  that  if  the  system  in  (A-l) 
is  controllable ,  then  the  transformation  matrix,  K,  is 
determined  by 


~-l 
K  =  E  E  x  , 


(A-8) 


where 
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;-i 


n 


3 

aH      .... 

n 

1 

4 

a5       .... 

1 

0 

5 

a6       .  .  .  . 

0 

0 

• 

• 

• 

• 

• 

• 

• 

• 

• 

1 

0       

0 

0 

0 

0       .... 

0 

0 

(A-9) 


and 


a  = 


n-1 

a 
n 


=  -E"1  An  B 


(A-10) 


The  elements  of  a  are  the  coefficients  of  the  characteristic 
polynomial : 


n      .  , 

det[SI-A]  =  det[SI-A]  =  S   +  Z      a.S1"1 

l=1   i 


(A-ll) 


The  matrix  inversion  in  equation  (A-10)  can  be  avoided  by 
using  the  Leverrier-Fadeev  method  for  calculating  the 
coefficients  of  the  characteristic  polynomial.   Once  the 
coefficients  are  known,  E    is  written  by  inspection. 
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Rane  [24]  presented  a  simplified  procedure  for  finding 
the  transformation  matrix,  K,  requiring  no  matrix  inversions. 
Substituting  equation  (A-4)  into  (A-l)  and  premultiplying 
equation  (A-2)  results  in: 


x  =  AKv  +  Bu 


=  KAv  +  KBu 


(A-12) 
(A-13) 


Comparison  of  equations  (A-12)  and  (A-13)  yields 


AK  =  KA 


(A- 14) 


and 


B  =  KB   . 


(A-15) 


Partition  K  into  n  column  vectors,  each  (nxl),  so  that 


K  —  L  k  -,   x0  ••••  k  J  • 

~1   ~  I  ~n 


(A-16) 


Substitution  of  (A-3)  and  (A-16)  into  (A-14)  and  (A-15)  gives 


^  ^1   ^2   -o       -n 


CSi  *2  !53 


k  ] 

~n 


1  .  .  .  .   0 
0  0 


0     0 


"al   -a2"--  an-l  'an 


(A-17) 
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and 


B  =  [k-  k„  k0  k  ] 

~1  ~2  ~3       ~n 


=  k   . 

~n 


(A-18) 


From  (A-17)  and  (A-18) 


k   =  B 
~n 


k   ,  =  A  k   +  k   a 
~n-l   ~  ~n   ~n  n 


k   0=Ak   n+k   a   n 
~n-2    ~  ~n-l    ~n   n-1 


K. n     —    r\    K.  ~  K.        a« 

~ Z         ~  ~ 3    ~n   o 


kn  =  A  k0  +  k   a0  , 
~1    ~  ~ I         ~n   2 


(A-19) 


or,  in  general, 


k   ,  =  A  k   •  ,-,  +  k   a   .  ,, 
~n-i    ~  ~n-i+l    ~n   n-i  +  1 


(A-20) 


for  ix[l,  2,  ...,  n-1].   The  column  vectors  k, ,  ...,  k   are 

~1      '  ~n 

found  in  a  simple  recursive  manner  and  completely 'determine 
the  transformation  matrix. 
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APPENDIX  B 
INVERSION  OF  CAUER  I  AND  CAUER  II  FORMS 


This  program  was  written  in  FORTRAN  IV,  and  requires 
minimal  input.   The  only  information  required  is: 

1.  the  order  of  the  system 

2.  which  inversion  is  required 

3.  the  quotients  from  the  continued  fraction  expansion 
Multiple  data  sets  are  possible,  and  input  in  the  following 
format : 

Card   Columns  Description  Format 


1-3 


M  =  the  desired  order  transfer 
function 


13 


1-20 

21-40 
41-60 
61-80 


h, ,  the  first  quotient  of  either 
Cauer  I  on    Cauer  II  continued 
fraction  expansion 

h«,  the  second  quotient 

h~ ,  the  third  quotient 

h,  ,  the  fourth  quotient 


D20.13 

D20.13 
D20.13 
D20.13 


N 


1-20 

21-40 

41-60 
61-80 


hUM   7,    the    (4N-7)th   quotient 
such   that    4N-7j<2M 

\N_6>  W-6<2M 
h4N-5'  4N"5I2M 
\N-4'    4N-4i2M 


D20.13 

D20.13 
D20.13 
D20.13 


Four  quotients  per  card  until 
2*M  quotients  have  been  input, 
where  M  is  the  system  order. 
Assume  this  is  the  Lth  card. 
The  (L+l)th  data  card  begins 
the  second  data  set. 
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L+l     1-3      M=system  order  13 

4-6      K=l  for  Cauer  I  and  K=2  for  13 

Cauer  II  inversion 


The  computer  program  has  been  written  to  handle  up  to 
20th  order  systems  (M_<20).   If  it  is  required  to  work  with 
higher  order  systems,  only  one  card  change  must  be  made. 
The  specification  statement  is  modified  to  read: 

REAL* 8  A(N),  B(N),  C(N)/N*0./.,  D(N)/N*0./,  DZERO 

where  N  is  an  integer  no  larger  than  999.   This  restriction 
can  be  lifted  by  changing  statement  2  to  read: 

2  FORMAT (2IR) 

where  R  is  the  mantissa  of  log, n(N).   The  REAL*8  in  the 
specification  statement  indicates  that  all  following  varia- 
bles and  arrays  are  real  valued  and  in  double  precision. 
Modification  to  either  single  or  extended  precision  would 
require  changes  in  all  format  statements.   If  this  is 
desireable,  the  user  should  consult  references  [2  5]  and  [26]. 

Execution  time  has  shown  to  be  less  than  .18  seconds 
for  systems  of  order  10  or  less. 
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//  EXEC  FORTCLG 
//FORT.SYSIN  OD  * 

REAL*8  A( 20) ,B(20) , C ( 2  01/20*0 ./ » D (20 ) /20*0 . / , DZERO 
C 

C      ..•  REAO  IN  SYSTEM  ORDER,  £  WHICH  CAUER  INVERSION  , 
C 

16  READ(5,1)  M,K 

1  F0RMAT(2I3> 
C 

C      ...  REAO  IN  QUOTIENTS  FROM  C3NT.  FRAC.  EXPANSION  ., 
C 

READ<5,2)  <A( I),B( I),I=l,M) 

2  FORMAT(4D20.13) 
C 

C      ...  DETERMINE  IF  INVERSION  IS  CAUER  I  OR  CAUER  II  , 
C 

IFU.EQ.l  )  GO  T3  10 
C 

C  ...    CAUER    II     INVERSION    ... 

C  ...     INITIALIZATION    ... 

C 

C(M)    =    B(M) 

Ml    =    M-l 

D(M1)    =    1.0 

D(M)    ~    A(M)*C(M) 
C 

C  ...    ITERATION    ... 

C 

DO    3    1=1, Ml 

L    =    1  +  1 

K   =    N-L 

DO    4   J  =  1,L 

KJ    =    K+J 

KL    =    KJ-1 

KP    =    KJ+l 

C(KJ)    =    3(K)*D(<J) 

IF(J.NE.L)    C(KJ)    =    C(KJ)+C(KP) 

IF(KL.EQ.O)     GO   TO    5 

O(KL)    =    A(K)*C<<L)+D(KJI 

GO    TO   4 
5    DZERO   =    1.0 
4    CONTINUE 

D(M)    =    A(K+1)*C<M) 

3  CONTINUE 
GO  TO  11 

C 

C      ...  CAUER  II  INVERSION  ... 

C      ...  INITI  ALIZATI3M  ... 

C 

10  C(M)  =  BIM) 

MM  =  M-l 

D(MM)  =  MM)*CM) 

D(M)  =  1.0 
C 
C      ...  ITERATION  ... 


C 


DO  6  1  =  1,  MM 

IP  =  1+1 

MI  =  M-IP 

DO  7  J  =  l,  IP 

MJ  =  MI+J 

ML  =  MJ-1 

MP  =  MJ  +  1 

CtMJJ  =  3(MI)*D(MJ)+C(MJ) 

IFCML.EQ. 0)  GO  TO  8 

D(ML)  =  A(MI)*CMJ)+DML) 

GO  TO  7 
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8  OZERO  =  1.0 

7  CONTINUE 

6  CONTINUE 
C 

C      ...  WRITE  OUT  TRANSFER  FUNCTION  IN  TERMS  OF 
C  NUMERATOR  AN3  DENOMINATOR  COEFFICIENTS 

C  WITH  APPROPRIATE  P3WER  OF  S  ... 

C 

11  WRITE(6,12) 

12  FORMAT (///l IX, «NJMERATOR« , 15X ,  •  DENOMINATOR  •  , 
110X, 'POWER  OF  S» ) 

WRITE(6,13)  DZER3,M 

13  FORMAT(///31X,D20.13,IOX, 12) 
DO  14  1=1, M 

MM I  =  M-I 

WRITE(6,15)  C(I),D( I),MMI 

14  CONTINUE 

15  FORM AT <// /6X,D2D.l 3, 4X, 020.13,1 OX, 12) 
C 

GO  TO  16 
END 
C 

C      FOR  ACTUAL  RUN  THIS  CARD  IS   /*   IN  COLUMNS  1  AND  2 
C 

//GO.SYSIN  DD  * 
C 

C      DATA  INPUT 
C 
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APPENDIX  C 

DETERMINATION  OF  WEIGHTING  MATRIX 
(Q)  FOR  PRESCRIBED  EIGENVALUES 

This  FORTRAN  IV  program  was  used  exclusively  on  the 
Naval  Postgraduate  School's  IBM  360/67  digital  computer 
and  includes  the  associated  job  control  language  statements 
The  program  consists  of  a  main  program  and  nine  subroutine 
subprograms.   The  purpose  of  each  subroutine  is  delineated 
below. 


SUBROUTINE 


DESCRIPTION 


READ 

RAMAT 

MULTPH 

HMATRX 

POLYNM 
HELP 

QTILDA 

QPIJ 

WRITE 


read  in  coefficients  of  numerator  and  denomi- 
nator polynomials  of  both  transfer  functions, 
and  places  each  system  in  phase  variable  form. 

determines  the  Routh  array  matrix  and  the 
transformation  matrix,  P. 

multiples  two  matrices,  Y  and  Z,  and  gives 
the  resulting  matrix  YZ.~ 

determines  the  quotients  of  continued  fraction 
expansion,  H1(H2),  and  the  state  matrix  in 
Cauer  II  f orm,~HHl(HH2 ) . 

T 
determines  the  product  det (SI-A)xdet (SI+A:.  ) 

computes  the  matrices  G,  T,  and  W  as  given 
in  Chapter  IV. 

computes  the  matrices,  Q,  and  Q  from  results 
of  subroutine  HELP. 

determines  the^off  diagonal  elements,  q. . , 
of  the  matrix  Q.  -1 

writes  all  two-dimensional  matrices. 
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The  input  required  has  been  reduced  to  a  minimum.   Multiple 
data  sets  are  possible;  and  input  as  indicated  below: 

Card   Columns  Description  Format 

1  1-3     N;  the  system  order  13 

2  1-16,     a   _  for  ie[l , 2 , . . . ,n] ;  F16.8 
17-32,     tne  denominator  coefficients  of 

33-48     the  known  transfer  function, 
for  an  Nth  order  system,  L 
cards  are  required,  where  L 
=K+1  and  K  is  the  integral  part 
of  N/5. 

L+2    1-16,     bn_1  for  i.e[l ,  2  ,  .  .  .  ,n]  ;  the  F16  .  8 

17-32,     numerator  coefficients  of  the 
3 3-48     known  transfer  function,  L 
cards  required. 

2L+2    1-16,  a    for  is [1 , 2 , . . . ,n] ;  F16.8 

17-3  2,  tne  denominator  coefficients 

33-48,  of  the  transfer  function  with 

...  prescribed  eigenvalues. 

3L+2    1-16,  3     for  is [1 , 2 , . . . ,n] ;  the  F16 . 8 

17-32,  numerator  coefficients  of  the 

33-48,  transfer  function  with  pre- 

...  scribed  eigenvalues. 

for  multiple  data  sets,  repeat 
the  same  prodecure.  Each  data 
set  requires  4L+1  cards. 


This  program  has  been  written  to  accept  systems  up  through 
20th  order.   To  increase  the  capability  of  the  program,  only 
the  dimension  statements  and  the  second  continuation  card  of 
the  equivalence  statement  require  modification.   The  system 
order  capability  can  be  increased  to  50.   Beyond  50th  order, 
the  program  requires  an  excessive  amount  of  storage  space 
(>510K  bytes).   Even  this  limitation  is  easily  overcome  by 
removing  the  four  cards  between  statements  1000  and  1001. 
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Any  other  modifications  (i.e.,  single  or  extended  pre- 
cision) require  extensive  changes  to  all  subprograms.   In 
this  case,  the  user  should  consult  [25]  and/or  [26].   It  is 
recommended  that  an  object  deck  or  disk  storage  be  used  when 
available  as  compilation  time  is  approximately  70-80%  of 
total  CPU  time. 
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